AD-A058  413  AIR  FORCE  INST  OF  TECH  WRIGHT-PATTERSON  AFB  OHIO  F/G  12/1 

AN  ADAPTIVE  FINITE-DIFFERENCE  METHOD  FOR  SOLVING  THE  STURM-LIOU— ETC (U) 
MAY  78  D A NELSON 

UNCLASSIFIED  AFIT-CI-78-22  ML 


0684  13 


7-  Cl-  W-Jx 


The  Pennsylvania  State  University 
The  Graduate  School 
Department  of  Computer  Science 


An  Adaptive  Finite-Difference  Method  for 
Solving  the  Sturm-Liouville  Problem  * 


A Thesis  in 


Computer  Science 


David  Alan 


Submitted  in  Partial  Fulfillment 
of  the  Requirements 
for  the  Degree  of 


Doctor  of  Philosophy 


born 


document 

>ub’ic  re*'- 

i out  ion  is  unlimited 


j 


SECURIT\PCi.  ASSIGNATION  OF  THIS  PAGE  fWi an  Data  Enlared; 


READ  INSTRUCTIONS 
BEFORE  COMPLETING  FORM 


. REPORT  NUMBER  12.  GOVT  ACCESSION  NO.I  3.  RECIPIENT'S  CATALOG  NUMBER 


REPORT  DOCUMENTATION  PAGE 


. REPORT  NUMBER 

Cl  78-72 


«.  TITLE  rand  Subtitle) 

An  Adaptive  Finite-Difference  Method  of 
Solving  the  Sturm-Liouville  Problem 


7.  AUTHORS 

Captain  David  A.  Nelson 


9.  PERFORMING  ORGANIZATION  NAME  AND  ADDRESS 


S.  TYPE  OF  REPORT  « PERIOO  COVERED 


Dissertation 


6.  PERFORMING  ORG.  REPORT  NUMBER 


8.  CONTRACT  OR  GRANT  NUMBERfe; 


10.  PROGRAM  ELEMENT,  PROJECT.  TASK 
AREA  A WORK  UNIT  NUMBERS 


■ 

I 


AFIT  Student  at  the  Pennsylvania  State 
University,  University  Park,  PA 


II.  CONTROLLING  OFFICE  NAME  AND  ADDRESS  12-  REPORT  DATE 

1978 

AFIT/CI  13.  NUMBER  OF  PAGES 

WPAFB  OH  45433  178  Pages 


MONITORING  AGENCY  NAME  A ADORESSfff  dlllormt  from  Controlling  Office;  15.  SECURITY  CLASS,  f of  ffifa  report; 

Unclassified 


15a.  DECLASSIFICATION/ DOWNGRADING 
SCHEDULE 


16.  DISTRIBUTION  STATEMENT  (ol  ffife  Report) 


Approved  for  Public  Release;  Distribution  Unlimited 


17.  DISTRIBUTION  STATEMENT  (of  the  abetrect  entered  In  Block  20,  It  different  from  Report) 


19.  KEY  WORDS  (Continue  on  reveree  tide  It  neceeaery  end  Identify  by  block  number ) 


20.  ABSTRACT  ( Continue  on  reveree  aide  It  neceeamry  end  Identify  by  block  number) 


The  signatories  below  indicate  that  they  have  read 
and  approved  the  thesis  of  David  Alan  Nelson. 


Date  of  Signature i 


Charlotte  F.  Fischer 
Frofessor  of  Computer  Science 
Chairman  of  Committee 
Thesis  Advisor 


Patrick  C.  Fischer 
Head  of  the  Department  of 
Computer  Science 


Don  E.  Heller 
Assistant  Professor  of 
Computer  Science,  Member  of 
Doctoral  Committee 


Wendell  H.  Millsf,  Jr. 
Assistant  Professor  of 
Mathematics,  Member  of 
Doctoral  Committee 


Jeffrey  R.  Spirn 
Assistant  Professor  of 
Computer  Science,  Member  of 
Doctoral  Committee 


We  evaluate  the  performance  of  our  algorithm  in 
terms  of  space  and  time  requirements  versus  the  eigen- 
value index  k.  to  ensure  that  the  method  has  practical 


iv 


application.  We  find  that  space  and  time  do  not  increase 
prohibitively  as  k increases.  Finally,  our  algorithm  is 
compared  to  a current  initial  value  method.  We  illustrate 
the  greater  accuracy  and  speed  of  our  method. 

Numerical  solutions  are  presented  for  several 
problems,  both  regular  and  singular,  along  with  other 
numerical  results  used  to  substantiate  our  theoretical 
results. 


TABLE  OP  CONTENTS 

Page 


LIST  OF  TABLES ix 

LIST  OF  FIGURES x 

ACKNOWLEDGEMENTS  xi 

Chapter 

1.  INTRODUCTION  1 

2.  STURM-LIOUVILLE  THEORY  4 

2.1  Introduction  .....  4 

2.2  Regular  Problems  4 

n 

2.3  Singular  Problems  8 

3.  EXISTING  METHODS  12 

3.1  Introduction 12 

3.2  Initial  Value  Methods 12 

3.3  Function  Space  Approximation  Methods  ....  14 

3.4  Conversion  to  an  Integral 17 

3«5  Separation  of  Variables  .....  18 


n 

11 


3.6  Approximation  of  Coefficients  21 

4.  FINITE-DIFFERENCE  METHODS  24 

4.1  The  Stormer  Uniform  Mesh  Method  ......  24 

4.1.1  Truncation  Error  26 

4.1.2  Convergence 26 

4.1.3  Solution  of  the  Matrix  Problem 27 

4.1.4  Error  Estimate  and  Deferred  Correction  . 29 


TABLE  OP  CONTENTS 


Chapter  Page 

4.2  The  Numerov  Uni  form  Mesh  Method 31 

4.2.1  Truncation  Error  ............  34 

4.2.2  Convergence 35 

4.2*3  Solution  of  the  Matrix  Problem  .....  36 

4.2.4  Error  Estimate  and  Deferred  Correction  . 37 

4.2.5  Numerical  Results  39 


4.3  A Non-Uniform  Finite-Difference  Method  . . . 4l 

4.3.1  Properties  of  the  Matrices  A and  B . . . 44 


4. 3. 1.1  Quasi -symmetry  of  A 4? 

4. 3. 1.2  Properties  of  B 49 

4.3.2  Truncation  Error  53 

4.3.3  Properties  of  the  Eigenvalues 54 

4.3.4  Convergence  .....  60 

4.3.5  Error  Estimate  ......  65 

4.3.6  Solution  of  the  Matrix  Problem  .....  69 

4.3.7  Deferred  Correction  77 

5.  AUTOMATIC  METHODS  80 

5*1  Introduction 80 

5.2  A Non-adaptive,  Iterative  Method  ......  81 

5.2.1  Background 81 

5.2.2  Truncation  Error 85 

5.2.3  The  Non-adaptive,  Iterative  Algorithm  . 91 

5.2.4  Numerical  Results  93 


TABLE  OF  CONTENTS 


vii 


Chapter  Page 

5.3  An  Adaptive,  Iterative  Method  96 

5.3.1  Refinement  Strategy  98 

5.3.2  Refinement  Criteria  102 

5.3. 2.1  Equidistribution  Strategy  102 

5. 3. 2. 2 Solution-Weighted  Strategy  103 

5. 3. 2. 3 D“*  Solution-Weighted  Strategy  . . . 104 

5.3.3  The  Adaptive.  Iterative  Algorithm  . . . 106 

5.3.^  Numerical  Results  109 

6.  STARTING  VALUE  AND  RECYCLING  PROBLEMS  117 

6.1  Starting  Value  Problem  117 

6.1.1  Objectives  of  the  Algorithm  ......  117 

6.1.2  Sturm  Sequences 119 

6.1.3  Bisection  Algorithm  124 

6.1.4  Starting  Value  Algorithm  127 

6.1.5  Numerical  Results  128 

6.2  Recycling  Problem  130 

7.  SINGULAR  PROBLEMS  135 

7.1  Introduction 135 

7.2  Infinite  Domain 136 

7*3  Singularity  in  q(x) 140 

7.4  Numerical  Results  142 

7.5  Order  of  Convergence 1 51 

7.6  Calculation  of  a Transformation  Function  . . 152 


viii 


TABLE  OF  CONTENTS 

Chapter  Page 

8.  PERFORMANCE  EVALUATION 154 

8.1  Space  Requirements . 154 

8.1.1  Regular  Problems  154 

8.1.2  Singular  Problems  157 

8.1.3  Conclusions  of  Space  Requirement 

Findings 158 

8.2  Time  Requirements 1 58 

8.3  Comparison  With  an  Initial  Value  Method  . . 170 

9.  SUMMARY  AND  CONCLUSIONS 173 

REFERENCES I76 


ix 


LIST  OF  TABLES 


Table  Page 

4.1  Results  for  kQ 40 

4.2  Results  for  k^ 41 

4.3  Convergence  Results  for  71 

4.4  Convergence  Results  for  k^ 71 


5.1  Results  for  Problem  I,  kQ  With  Algorithm  5*4  . . 94 

5.2  Results  for  Problem  I,  k^  With  Algorithm  5.4  . . 94 

5.3  Results  for  Problem  II,  kQ  With  Algorithm  5.4  . 95 

5.4  Results  for  Problem  II,  1 2 With  Algorithm  5.4  . 95 


5«5  Comparison  of  Refinement  Strategies  106 

5.6  Results  of  Algorithm  5*6  . Ill 

6.1  Results  of  Algorithm  6.3  on  Problem  I 129 

6.2  Results  of  Algorithm  6.3  on  Problem  III  ....  129 

6.3  Recycling  of  Ai 1 31 

6.4  Recycling  With  Algorithm  6.5 134 

7.1  Singular  Problems  144 

7.2  f'(x)Ax  for  Problem  C 153 

8.1  Adaptive  vs  SLEIGN  - Regular  Problems  171 

8.2  Adaptive  vs  SLEIGN  - Singular  Problems  Computed 

by  Both  Algorithms 171 

8.3  Adaptive  vs  SLEIGN  - All  Singular  Problems  . . . 172 


X 


LIST  OF  FIGURES 

Figure  Page 

4.1  The  Matrix  A 33 

5.1  Mesh  Before  Refinement 98 

5.2  Mesh  After  Refinement 98 

5.3  Point  x.  Before  Splitting  100 

J 

5.4  Point  x.  After  Splitting 100 

J 

5.5  Mesh  Distribution  for  Problem  III,  1 1 5 

5.6  Mesh  Distribution  for  Problem  III,  116 

?.l  Mesh  Distribution  for  Problem  B,  14? 

7.2  Mesh  Distribution  for  Problem  B,  148 

7.3  Mesh  Distribution  for  Problem  C,  149 

7.4  Mesh  Distribution  for  Problem  C,  * 150 

8.1  N vs  i - Problem  I 161 

8.2  N vs  i - Problem  II 162 

8.3  N vs  i - Problem  III 163 

8.4  N vs  i - Problem  B 164 

8.5  T vs  i - Problem  I I65 

8.6  T vs  i - Problem  II 166 

8.7  T vs  i - Problem  III 167 

8.8  T vs  i - Problem  3 168 

8.9  T vs  M - Problem  I . . I69 


I 


xi 

ACKNOWLEDGMENTS 

The  author  expresses  his  gratitude  to  Dr.  Charlotte 
F.  Fischer,  his  thesis  advisor,  for  education,  motivation, 
and  guidance  during  his  work  at  The  Pennsylvania  State 
University. 

Support  for  the  author's  graduate  work  has  been 
provided  by  the  United  States  Air  Force  through  the  Air 
Force  Institute  of  Technology. 


Chapter  1 


INTRODUCTION 

The  Sturm-Liouville  boundary  value  problem  arises 
frequently  in  the  study  of  mathematical  physics.  For 
example  the  problem 

y ♦ Ay  * 0 (1.1.1) 

y(0)  = y(l)  * 0 

occurs  in  the  study  of  the  vibrating  elastic  string.  The 
values  of  A for  which  a solution  exists  are  found  to  be 
related  to  the  natural  frequencies  of  the  string.  The 
problem  (1.1.1)  is  an  example  of  a regular  Sturm-Liouville 
problem. 

A singular  problem  such  as 

y ♦ C*  - qU))y  a 0 

y(0)  ■ 0,  y(x)  bounded  as  x *•  «° 


occurs  in  quantum  mechanics.  Examples  include  the  study 
of  nuclear  shell  theory  and  the  vibrational  theory  of  di- 
atomic molecules.  Regular  and  singular  problems  are  for- 


mally defined  in  Chapter  2. 

Many  methods  have  been  proposed  to  solve  the  Sturm- 


Liouville  problem  numerically.  Among  these  is  the  general 
class  of  finite -difference  or  direct  methods.  This  thesis 


2 


reviews  the  various  methods,  placing  particular  emphasis 
on  the  finite-difference  type.  An  automatic,  adaptive 
finite-difference  method  is  developed  in  which  the  mesh  of 
grid  points  used  in  the  discretization  depends  on  the 
solution. 

In  Chapter  2 we  provide  the  reader  with  necessary 
background  material  by  presenting  many  of  the  classical 
results  of  Sturm-Liouville  theory. 

In  Chapter  3 we  review  other  methods  for  solving 
the  Sturm-Liouville  problem  numerically.  The  theoretical 
basis  of  each  method  is  presented  along  with  discussions  of 
recent  work  in  each  area. 

The  review  continues  into  Chapter  4 with  Sections 
4.1  and  4.2  exploring  in  detail  the  uniform  mesh  finite- 
difference  nethod.  In  Section  4.3  we  develop  a non-uniform 
mesh  finite-difference  method.  Properties  of  the  resulting 
matrix  problem  are  established,  a solution  algorithm  is 
developed,  and  convergence  results  as  well  as  an  error 
estimate  are  derived. 

Chapter  5 deals  with  automatic  solution  methods. 

Both  a non-adaptive  and  an  adaptive  algorithm  are  presented 
along  with  numerical  results  from  each. 

For  a numerical  algorithm  to  be  completely  automatic 
the  user  can  be  expected  to  supply  only  the  problem  to  be 
solved.  He  must  not  be  required  to  provide  an  estimate  to 
the  solution.  To  this  end,  in  Chapter  6 we  present  an 
algorithm  for  computing  an  estimate  to  the  solution  for  use 


f.' 


as  a starting  value  in  the  algorithms  of  Chapter  5. 
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In  Chapter  7 we  develop  some  additional  techniques 
necessary  to  solve  singular  problems  adaptively,  ^wo  gen- 
eral types  of  singular  Sturm-Liouville  problems  are  solved 
and  numerical  results  presented.  For  each  problem  solved, 
the  adaptive  method  is.  compared  to  a uniform  mesh  finite- 
difference  method.  We  also  note  a further  use  of  the 
adaptive  method  showing  how  it  can  be  used  to  select  the 
best  transformation  function  for  use  with  some  uniform 
mesh  methods. 

Finally,  in  Chapter  8 we  evaluate  the  performance  of 
the  method  by  showing  the  relationship  between  the  index 
of  the  desired  eigenvalue  and  the  space  and  time  required 
to  compute  that  eigenvalue.  In  addition  we  compare  its 
performance  to  an  existing  software  package  developed  at 
Sandia  Laboratories,  Albuquerque,  New  MexicoLll. 


r 


Chapter  2 


STURM-LIOUVILLE  THEORY 


2.1  Introduction 

As  was  stated  in  the  previous  chapter  the  purpose 
of  this  thesis  is  to  study  adaptive  methods  for  solving 
numerically  a specific  type  of  Sturm-Liouville  system. 

To  solve  any  problem  numerically  it  is  first  nec- 
essary to  understand  the  mathematical  theory  of  the  prob- 
lem. To  this  end  we  present  in  this  chapter  an  overview  of 
the  theory  of  Sturm-Liouville  systems,  both  regular  and 
sineoilar . 

2 . 2 Regular  Problems 

A Sturm-Liouville  (S-L)  equation  is  a second-order, 
homogeneous,  linear  differential  equation  of  the  form 

dx  aJ]  L*r(x)  - q(x)]y  - 0,  a < x < b (2.2.1) 

where  ^ is  a parameter,  p,  r,  and  q are  real -valued  func- 
tions of  x.and  p and  r are  positive.  We  further  require 
that  p,  p , r,  and  q be  continuous  on  la,b]. 

A Sturm-Liouville  system  is  a Sturm-Liouville 
equation  together  with  boundary  conditions. 

A regular  S-L  system  is  an  S-L  equation  defined  on 
a finite  interval, La, b ],  with  boundary  conditions  of  the 


L 


r 


o 


i 


form 


5 


A1y(a)  * A2y  (a)  0 (2.2.2) 

BlV(b)  ♦ B2y'(b)  = 0 (2.2.3) 

where  A^t  A2,  B^,  and  B2  are  real  constants  and 

I Ax | + |A2I  / 0 
|B1I  ♦ |B2I  / 0 

The  values  of  V for  which  (2.2.1)  has  a non-trivial 
solution  satisfying  (2.2.2)  and  (2.2.3)  are  called  eigen- 
values . Any  solution  corresponding  to  an  eigenvalue  is 
called  an  eigenfunction . The  set  of  all  eigenvalues  is 
called  the  spectrum  of  the  system. 

We  note  here  that  any  equation  of  the  form  (2.2.1) 
can  be  transformed  via  the  Liouville  transformation  Lsee 
Birkhoff  and  Rotal_3].  p.  104,  for  example]  into  liouville 
normal  form 

- q(x)  ]w  0 

dx^ 

which  is  defined  on  the  new  interval  |_c,d].  We  therefore 
choose  to  work  primarily  with  S-L  equations  of  the  form 

y"  - O - q(x)]y  - 0.  (2.2.4) 

In  addition  we  choose  to  consider  only  boundary  conditions 
of  the  first  kind,  i.e.  conditions  (2.2.2)  and  (2.2.3) 
with  A 2 = B 2 — 0.  This  serves  to  simplify  the  matrix 


r 
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0 

n 


eigenvalue  problems  we  must  solve i however,  extensions  to 
other  boundary  conditions  are  possible.  FoxLlBj  discusses 
some  of  these  techniques  as  they  relate  to  two-point  bound- 
ary value  problems. 

We  are  now  able  to  state  the  following  results 
about  S-L  systems.  These  results  are  presented  here  as 
theorems,  however  all  proofs  are  omitted.  The  reader  is 
referred  to  Hillel.24]  or  Coddington  and  Levinsonl.9]  for 
the  proofs. 

The  first  few  theorems  establish  the  existence  of 
solutions  and  some  basic  properties  of  these  solutions. 

THEOREM  2.1 » The  regular  S-L  system  (2.2.1), 

(2.2.2),  (2.2.3)  has  an  infinite  number  of  real  eigenvalues 

^0»  . *2’  f°rmin£  a monotone  increasing  sequence 

with  » as  n - 

n 

THEOREM  2.2«  The  eigenfunction  yn(x)  corresponding 
to  *n  has  exactly  n zeros  in  the  interval  La.b]  and  is 
uniquely  determined  up  to  a constant  factor. 

THEOREM  2.3»  Eigenfunctions  of  a regular  S-L  system 
corresponding  to  different  eigenvalues  are  orthogonal  with 
weight  function  r,  i.e. 

/ r(x)yB<x)yn(x)dx  * 0 IT  »m  / 

(Note:  In  equation  (2.2.4)  r(x)  *•  1.) 

The  next  several  results  snow  how  the  eigenfunctions 


l 
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corresponding  to  different  eigenvalues  differ  in  terras 
of  zeros  and  maxima-rainiraa. 

THEOREM  2.4«  Consider 

u (x)  + F(x)u(x)  - 0 

where  F(x)  is  positive  and  continuous  in  (a,b).  Let  u^(x) 
and  u2(x)  be  2 linearly  independent  solutions  of  (2.2.4) 
and  suppose  that  u^(x)  has  at  least  2 zeros  in  (a,b).  If 
x^  and  x2  are  2 consecutive  zeros  of  u^(x)  then  u2(x)  has 
exactly  one  zero  in  (x^,x2).  (i.e.  zeros  are  interlaced) 

THEOREM  2. 5»  Consider 

(a)  u (x)  ♦ F(x)u(x)  = 0 

(b)  v (x)  ♦ G(x)v(x)  = 0 

where  F(x)  and  G(x)  are  positive,  continuous  and  G(x)  > F(x) 
in  (a,b).  Suppose  that  (a)  has  a solution  with  2 consec- 
utive zeros  at  x ~ x^  and  x - x2  and  that  a < x^  < x2  < b. 
Then  if  v(x)  is  a solution  of  (b)  with  a zero  at  x - x^, 
v(x)  has  at  least  one  zero,  x^,  with  x^  < x^  < x2.  (i.e. 

the  distance  between  consecutive  zeros  decreases  as  X 
increases) 

THEOREM  2.6 i With  the  same  assumptions  on  ^(x) 

and  G(x),  let  u(x)  and  v(x)  be  solutions  of  (a)  and  (b) 

• • 

such  that  u(x^)  - v(Xj)  - 0,  u (x^)  v (x^)  ->  0.  If 
u(x)  is  increasing  in  (.x^tx2)  and  reaches  a maximum  at 


x = then  v(x)  reaches  a maximum  at  some  x_  with 

2 3 

xi  < x3  < x2. 

The  following  theorem  gives  us  a means  for  bound- 
ing the  ith  eigenvalue. 

THEOREM  2.7»  Consider  the  system  ( 2. 2.1) , ( 2. 2. 2) , 
(2.2.3)  and  let 

q = min  [q(x)«x  € [a,b]} 

Q = max  (q(x)«x  € La.b]} 

Let  the  system  with  q(x)  replaced  by  q have  eigenvalues 
and  with  q(x)  replaced  by  <4 , tvjj*  Then 

< vi  for  all  i* 

2.3  Singular  Problems 

A singular  S-L  system  is  an  S-L  equation  defined 
either  on  an  infinite  interval  or  on  a finite  interval 
which  has  a singular  point  of  the  function  q at  an  endpoint. 

This  definition  applies  to  the  normal  form  (2.2.4). 
When  the  general  equation  (2.2.1)  is  considered,  a sing- 
ular system  may  also  have  a singularity  or  zero  of  p or  r. 
The  results  which  follow  are  applicable  to  the  normal  form. 

We  consider  the  case  where  the  interval  is  (0,*>). 

As  noted  by  Titchmarsh(.4ol  the  case  of  a finite  interval 
with  a singularity  at  one  end  is  quite  similar. 
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THEOREM  2.8:  For  every  non-real  value  of  X the 

system  has  at  least  one  non-trivial  solution  in  L (0,°»)  . 

(i.e.  a solution  y(x)  such  that  ly  (x)dx  < *>) 

0 

This  type  of  S-L  system  can  be  placed  into  one  of 
two  classes.  If  all  solutions  are  in  L (0,°°)  we  have  the 
limit-circle  case  at  <».  If  only  one  solution  is  in  L2{0,°°) 
we  have  the  limit-point  case  at  «. 

THEOREM  2.9»  If,  for  a particular  value  of  X,  the 

p 

equation  (2.2.4)  has  2 independent  solutions  in  L (0,«) 
then  this  property  holds  for  all  X. 

An  important  sufficient  condition  for  the  limit- 
point  case  is  as  follows. 

THEOREM  2.10t  Suppose  q(x)  is  bounded  from  below. 
Then  the  limit-poinc  case  holds  at  <*>  and  the  solution  y(x) 
satisfies 

i)  lim  y(x)  y (x)  » 0 

X*'» 

ii)  y (x)  € L2(0,<») 

iii)  | q(x)  | ^y(x)  € L2(0,<*) 

As  with  regular  problems  we  can  say  something  about 
the  eigenfunctions. 

THEOREM  2.11»  Square -integrable  (i.e.  € L2(0t°<>)  ) 
eigenfunctions  u and  v belonging  to  different  eigenvalues 
of  a singular  S-L  system  are  orthogonal  with  weight  r 


1C 


1 

i] 

whenever 

L u(x)v'(x)  - v(x)u’(x)  J = 0. 

(3-b  x=a 

The  set  of  values  of  X for  which  the  limit-point 
case  holds  is  called  the  spectrum  of  the  S-L  system. 

THEOREM  2.12:  The  spectrum  is  a closed,  infinite, 
and  unbounded  set. 

The  spectrum  can  either  be  a pure  point  spectrum 
if  all  points  are  isolated,  a continuous  spectrum,  or  a 
combination  of  the  two. 

In  a regular  S-L  system  we  found  that  we  always  have 
a pure  point  spectrum,  but  in  a singular  system  the  nature 
of  the  spectrum  depends  on  the  function  q. 

THEOREM  2. lit  If  q(x)  * °°  with  x then  there  is  a 
pure  point  spectrum  with  at  most  one  spectral  value  to 
the  left  of  qQ  = inf  q(x). 

THEOREM  2. 14 1 If  q(x)  € L(0,°°)  then  there  is  a 
continuous  spectrum  covering  the  whole  positive  real  axis. 

There  may  be  a point  spectrum  on  the  negative  real  axis 
which  may  have  x = 0 as  a cluster  point. 
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the  continuous  spectrum  covers  (-  °°,0)  and  there  may  be 
a point  spectrum  on  (0,<«)  ; otherwise,  if  |q(x)|~^  / L(0,°°), 
the  continuous  spectrum  covers  the  whole  real  axis. 

Analysis  of  the  spectrum  of  a singular  problem  is 
usually  done  on  an  individual  basis  for  each  function  q(x). 

In  Chapter  7 we  discuss  singular  problems  further 
and  give  examples  of  each  type. 


o 

0 

n 
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Chapter  3 


EXISTING  METHODS 


3 • 1 Introduction 

The  task  of  solving  the  Sturm-Liouville  problem  has 
been  the  subject  of  much  research.  There  have  been  many 
numerical  methods  proposed  and  these  methods  are  vastly 
diverse  in  technique  and  theory.  We  can,  however,  cate- 
gorize the  methods  into  the  following  6 general  areas. 

a.  Initial  value  methods 

b.  Function  space  approximation 

c.  Conversion  to  an  integral 

d.  Separation  of  variables 

e.  Approximation  of  coefficients 

f.  Finite-difference  methods 

This  thesis  is  directed  toward  the  final  area,  that 
of  finite-difference  methods.  However,  in  order  to  provide 
the  proper  background  for  our  work  we  present  here  a brief 
description  of  areas  a.-e.  In  Chapter  4 we  present,  in 
greater  detail  than  is  demanded  here,  a description  of  the 
finite-difference  methods. 

3 . 2 Initial  Value  Methods 

Initial  value  or  shooting  methods  have  received  con- 
siderable attention  from  researchers,  especially  in  the 
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applied  areas  where  specific  eigenvalue  problems  are  being 
solved.  This  is  probably  due  to  the  fact  that  initial 
value  theory,  on  which  these  methods  rely  heavily,  is  wide- 
ly understood  and  good  code  for  solving  initial  value  prob- 
lems is  in  abundance. 

The  solution  of  the  S-L  system  (2.2.1)  and  (2.2.2) 
is  related  to  an  initial  value  problem  in  the  following  way. 
For  a fixed  value  of  1 we  consider 

Ly  = Ar(x)y  , a < x < b (3*2.1) 

A1y(a)  ♦ A2p(a)y  (a)  = 0 (3.2.2) 

Cj^yfa)  + C 2p ( a ) y (a)  = 0 

• • 

where  Ly  = -(p(x)y  ) q(x)y  and  where  and  C?  are  arbi- 
trarily chosen  constants  such  that  (A2C^  - C2A1^  ^ °* 

We  are  assured  that  (3*2.1),  (3*2.2)  has  a unique 
non-trivial  solution.  If  this  solution  is  denoted  by 
y(Ajb)  then  we  can  define  the  following  function. 

F(M  = B1y(XJb)  ♦ B2p(b)y  (Ajb) 

Requiring  that  F(A)  = 0 is  thus  equivalent  to  forc- 
ing our  solution  to  satisfy  the  boundary  condition  (2.2*3). 
So  we  see  that  F(A  ) = 0 if  and  only  if  *n  is  an  eigenvalue 
of  (3*2.1)  and  y(*nix)  is  the  corresponding  eigenfunction. 
The  S-L  problem  has  thus  been  reduced  to  finding  the  roots 
of  F(A)  - 0 and  the  corresponding  solutions  of  (3.2.1). 
Kellert25]  advocates  using  a fixed  point  iteration  scheme 
for  finding  the  roots.  Dranofft.l5j  solves  a mass  and  heat 


ry 
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transfer  problem  using  a shooting  method  and  Newton's 
method.  A more  slowly  converging  method  may  be  to  use  the 
bisection  method. 

There  are  many  examples  in  the  literature  of  each 
such  method  as  well  as  others.  Every  method,  however,  must 
solve  the  initial  value  problem  (3*2.1),  (3*2.2)  to  perform 
the  root  finding  procedure.  The  designer  of  a shooting 
method  algorithm  may  choose  any  of  the  numerous  initial 
value  methods  to  solve  the  problem.  Thus  we  see  that 
initial  value  methods  for  solving  the  Sturm-Liouville 
problem  are  varied  in  technique  and  are  dependent  both  on 
the  problem  to  be  solved  and  the  designer  of  the  method. 

3 . 3 Function  Space  Approximation  Methods 

The  class  of  methods  receiving  the  most  attention 
by  current  researchers  is  that  of  function  space  approx- 
imations. The  procedure  is  to  approximate  the  solution  of 
the  eigenvalue  problem  by  a linear  combination  of  linearly 
independent  functions  in  some  appropriately  chosen  function 
space.  The  coefficients  are  determined  by  requiring  that 
some  measure  of  the  error  be  minimized. 

For  all  methods  described  here  we  let  lu^(x)}<?_^  be 
the  set  of  basis  functions  for  our  space  of  admissible 
functions,  H. 

An  approximation  to  the  solution  of  the  eigenvalue 
problem  can  be  written  as 

N 

yN(x)  •-*  £ c.  u.  (x) 
i-0  1 1 


(3*3*1) 
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The  methods  to  be  discussed  here  differ  in  essen- 
tially two  ways.  First,  they  differ  in  the  way  in  which 
the  { c^ ) are  determined  and  second  in  the  choice  of  the 
space  H.  We  describe  two  of  the  more  common  ways  to  choose 
the  (Ci).  The  possible  choices  of  H are  evident  from  the 
cited  examples. 

i)  Rayleigh-Ritz  (Galerkin) 

Rayleigh  noted  that  the  eigenvalues  and  eigen- 
functions of  (2.2.1)  are  stationary  values  and  critical 
points  for 

R [u]  = N [u]  / D Lu] 

^ • 2 p 

where  N Lu]  = /(pu  + qu  )dx 

a 

b P 

and  D [u]  = S (ru  )dx 

a 

The  Rayleigh-Ritz  method  finds  the  critical  points 
and  values  of  D Lu]  on  some  finite-dimensional  subspace 
of  H,  where  H is  the  space  of  all  functions  u such  that 
D [u]  < °°. 

This  procedure  results  in  a matrix  problem 

Ac  = .ABc  (3-3.2) 

where 

b . . 

(a^j)  •=  .Mpu^Uj  * qUj^UjJdx 

St 

b 

lbij>  • 

A 

The  extrema  of  the  Rayleigh  quotient  R Lu]  are  at 


r. 
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the  yN(x)  - Ec^u^  wh-*se  coefficient  vectors  c are  the  eig- 
envectors of  (3.3.2). 

Birkhoff,  deBoor,  Swartz,  and  Wendrofft4]  use  both 
the  2n-dimensional  cubic  Hermite  subspace  and  the  (n+1)- 
dimensional  cubic  spline  subspace.  They  found  that  both 
methods  approximate  eigenvalues  with  error  O(h^)  and  eigen- 
functions  with  0(hJ), where  h = maxfx^  - with  (x^} 

being  the  partition  on  which  the  subspace  is  defined. 

Earlier  CourantCll]  suggested  using  piecewise  linear 
functions  as  the  subspace.  Birkhoff  and  Fix£5j  give  an  ex- 
ample using  the  trigonometric  polynomials  where  u^  = 1, 
u ^ = cos(nx),  and  u^.^  = sin(nx). 

In  each  case  the  matrix  entries  must  be  computed  by 
evaluating  integrals  and  thus,  while  the  method  can  be 
quite  accurate,  it  can  also  be  time  consuming. 

ii)  Collocation 

In  these  methods  the  approximate  solution  (3.3.I) 
is  required  to  satisfy  the  differential  equation  at  N 

|t 

distinct  points  including  the  boundary  conditions,  i.e. 

Lyjy (xi)  * Xr(x^) y^j (x^),  i=2i 3 » • • • |N-1 

y„(*l)  = yN<*„>  = 0 


The  resulting  matrix  eigenvalue  problem  to  solve 
for  the  coefficients  {c^}  is 

(A  - XB)c  = 0 
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where 


<ai j*  * Luj<xi) 
<bu>  * rUjlUjUj) 


Russell  and  Varah£38]  discuss  subspaces  of  (2n-l) 
degree  B-splines  and  (2n-l)  degree  Hermites.  Russell  and 
Shampine[37]  show  some  results  using  general  piecewise 
polynomials. 

3.4  Conversion  to  an  Integral 

Another  technique  for  the  S-L  problem  is  to  convert 
the  differential  equation  into  an  integral  equation.  This 
involves  finding  the  Green's  function.  g(x,v),  for  the 
system  (2.2.1),  (2.2.2).  The  function  g(x,v)  is  such  that 
the  solution  of  (2.2.1)  is  given  by 


b 

y(x)  = l/g(x, v )r(v )y(v )dv 
a 
b 

or  y (x)  = X/G(x,v)y(tf )dv . (3.4.1) 

a 

If  a mesh  ivj}»  jsl,2,...,J,  is  now  introduced  on 
La, b], with  a < < v2  < •••  < * j i b,  we  can  define  o 

quadrature  rule 


J 

QtCf(v)}  = I a,  f ( v . ) 

J i=l  1 1 

to  approximate 

b 

Q[f(v)}  = f f(v)dv  . 


If  this  quadrature  rule  is  used  to  approximate  y(x^) 
by  y^  at  each  point  we  get  the  following  matrix 


t- 


eigenvalue  problem. 
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J 

yi  =A.£  <*iG(xivj)y<j'  .« 

j ~ ) 

...J 

or 

*L~1y  = Ay 

where 

y = • • -Yj)* 

and 

Aij  = aiG(xi,wj)* 

The  problem  is  now  solved  as  a matrix  eigenvalue 
problem  where  the  matrix  A is  non-symmetric  and  is,  in 
general,  a full  matrix.  Keller  shows  that  for  every  A of 
the  integral  equation  (3*^.l)  there  exists  an  eigenvalue 
of  the  matrix  A,  correct  to  within  an  accuracy  pro- 
portional to  the  error  in  the  quadrature  formula  used. 

Many  excellent  adaptive  quadrature  rules  for  in- 
tegration are  available  (see  Davis  and  Rabinowitzl.12]  for 
examples).  If  such  a rule  were  used  to  approximate  Q{  f («')}, 
the  resulting  S-L  method  could  be  considered  adaptive. 

Use  of  this  method  entails  computing  J integrals 
and  J constants,  and  the  solving  of  a matrix  problem  in 
which  A is  a full  matrix.  For  most  problems  this  technique 

would  be  highly  time  consuming  because  the  work  is  pro- 
2 

portional  to  J . It  is  for  this  reason  that  we  look  to 
finite-difference  methods  for  an  adaptive  algorithm.  By 
selecting  the  proper  discretization  and  the  best  solution 
methods,  the  amount  of  work  can  be  made  proportional  to  J. 

3.5  Separation  of  Variables 

This  method  is  based  on  the  Prufer  transformation 


ifliMHii 
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which  transforms  the  original  S-L  equation  (2.2.1)  into 
two  first-order  equations,  one  of  which  is  independent  of 
the  solution  of  the  other.  Two  new  variables,  t(x)  and 
0(x),  are  introduced  as  follows » 

y(x)  = t(x)sin  9(x) 
p(x)y  (x)  » t(x)cos  © (x) 

Under  this  transformation  we  get,  from  the  original 
problem, 

d0(x)/dx  = cos20(x)/p(x)  ♦ [q(x)  ♦ Xr (x) ]sin20(x) 

(3-5.1) 

subject  to  the  boundary  conditions 

9(a)  = A (3.5.2) 

0(b)  = B + (n-l)rr  ( 3 • 5 - 3 ) 

0 < A < n,  0 c B < rr , na  non-negative  integer. 

and 

dt(x)/dx  ~ t(x)sinO(x)cos0(x)Lp^y  - q(x)  - Xr(x)l 

(3-5.4) 

If  0(x»X)  denotes  the  solution  of  (3.5.1),  satis- 
fying (3-5.2),  then  the  eigenvalues  of  the  original  problem 
are  those  values  of  X for  which 

0(b»X)  = B ♦ (n-l)n. 


This  means  that  the  calculation  of  the  nth  eigen- 
, * r is  equivalent  to  finding  the  solution  of 
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B 4-  (n-l)tr  - e(bjX)  = 0.  (3.5.5) 

If  we  define  X(x,X)  = d0(x;X)/dX,  then  differenti- 

ating (3.5.1)  we  get 

”™dxA^  = X(x»X)Lq(x)  * Xr(x)  - ]sin28 (x;* ) 

+ r(x)sin2e(x;M  (3.5*6) 

X(a,X)  = 0 (3.5.7) 

and  we  can  solve  (3*5.5)  with  the  Newton-Raphson  method 

xik+1)  = Xik)  * * (n-l)w  - 0(b,X(k)) lX"1(b,X(k)) 

n n n n 

(3.5.8) 

The  numerical  algorithm  for  finding  the  (n+1) st 
eigenvalue  then  proceeds  as  follows: 

1.  Let  be  an  approximation  to  X 

Repeat  the  following  steps  for  k~0,l,... 

2.  Solve  (3. 5.1) -(3* 5. 2)  and  ( 3 . 5 .6) - ( 3 . 5 .7)  by 
some  reliable  ODE  solver. 

3.  Compute  X^k+1)  by  (3-5.8) 

4.  If  |X*k>1)  - \*k*|  ,<  <■-  STOP 

n n 

Banks  and  KurowskiL2]  recommend  that  the  inte- 
gration be  carried  out  as  follows: 

1.  Integrate  (3. 5 . 1 ) -(3 . 5. 2)  from  a to  (a-b)/2  to 
get  0a(x,k ) . 

2.  Integrate  (3. 5.1) -(3. 5.2)  from  b to  (a*b)/2 
to  get  6b(x,X) . 


We  now  want 
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EU  ) = 

n a d n 


0b  (^T'kn)  = 0 


and  a similar  Newton-Raphson  technique  is  used  to  solve 
this  equation. 

This  is  the  technique  that  Bailey,  Gordon,  and 
5hampine[l]  use  in  their  recent  paper.  The  method  can 
also  be  considered  adaptive  because  a variable  step  inte- 
grator can  be  used  to  solve  for  0„  and  0,  . 

a b 

Since  the  eigenvalues  are  unaffected  by  the  Prufer 
transformation,  the  accuracy  of  these  methods  depends 


entirely  on  the  accuracies  of  the  integration  and  root- 
finding methods  used. 


3.6  Approximation  of  Coefficients 


Canosa  and  Gomes  de  Oliveiral.8]  first  used  this 
technique  to  solve  the  one-dimensional  Schrodinger  equation 

d2y/dx  - V(x)y  + Ey  - 0 
y(0)  = y (L)  * 0 

The  method  is  based  on  the  idea  that  if  the  poten- 


tial, V(x),  is  replaced  by  a step  function,  a greatly  simp- 


lified problem  results  and  can  be  solved  exactly. 

i ' ' 


V(x)  is  approximated  as  follows: 


n yields  a 


Now  each  interval  (x^_^,x^),  i=l,2 
problem  of  the  form 

d2y/dx2  <•  (E  - V.)y  = 0. 

It  is  well  known  that  such  an  equation  has  a basis 
of  solutions  consisting  of  piecewise  sine  and  cosine,  or 
exponential  functions.  By  enforcing  the  boundary  con- 
ditions  and  continuity  of  y and  y , the  problem  reduces 
to  that  of  solving  a homogeneous  linear  system. 

The  proponents  of  this  method  claim  an  important 
advantage  over  many  other  methods.  Because  the  basis 
functions  depend  transcendentally  on  X,  the  system  which 
results  has  an  infinite  number  of  roots.  In  most  methods 
for  approximating  eigenvalues  it  is  only  practical  to 
approximate  a small  number  of  values  with  any  given  mesh 
size.  For  higher  eigenvalues  the  relative  error  becomes 
prohibitively  large.  In  the  approximation  of  coefficients 
methods  the  relative  error  is  much  more  uniform  as  a func- 
tion of  the  eigenvalue  index,  k.  This  enables  Preuss[35l, 
for  example,  to  solve  one  problem  for  *10o  his  method 

to  a relative  accuracy  of  6.0x10"^.  Such  an  undertaking 
would  be  impractical  for  most  other  methods. 

Preuss  extends  the  above  work  in  two  directions.  He 
considers  the  general  S-L  problem  and  he  generalizes  Canosa 
and  Gomes  de  Oliveira's  work  to  include  mth  degree  poly- 
nomial approximations  of  the  coefficient  functions.  In- 
stead of  solving  the  resulting  simplified  problem  exactly 


he  uses  an  initial  value  method  based  on  Taylor's  series 
expansions,  thus  eliminating  the  need  for  knowing  the  basis 
functions  exactly. 

Canosa  and  Gomes  de  Oliveira  found  their  method  to 

be  0(h)  for  both  eigenvalues  and  eigenfunctions,  where 

h=  max(x^4^'-  x^).  Preuss  shows  that  for  mth  degree  ap- 
i m+ 1 

proximations  the  order  is  0(h  ). 

The  final  category  of  solution  methods  is  finite- 
difference  methods.  In  Chapter  4 we  present  these  methods 
and  introduce  a non-uniform  method  which  we  use  to  solve 
S-L  problems  adaptively. 


Chapter  4 


FINITE -DIFFERENCE  METHODS 

In  this  chapter  we  describe  the  finite-difference 
method  in  detail.  First,  in  Section  4.1  we  introduce  the 
simplest  uniform  mesh  method,  Stormer's  method,  and  ana- 
lyze it  completely.  Then,  in  Section  4.2,  we  examine  a 
higher  order  method  and  discuss  the  resulting  changes  in 
convergence  theorems,  solution  techniques,  and  error 
estimates.  In  Section  4.3  we  present  a non-uniform  mesh 
finite-difference  method  which  will  become  the  foundation 
for  our  adaptive  S-L  solver.  Results  pertaining  to  all 
the  important  properties  of  the  method  are  presented. 

4.1  The  Stormer  Uniform  Mesh  Method 

First  we  recall  the  problem  to  be  solved. 

y *•  - q(x))y  = 0 (4.1.1) 

y(a)  = y(b)  - 0 (4.1.2) 

where  q(x)  is  continuous  in  [a.bj.and  a and  b are  finite. 
Let  the  interval  [a,b]  be  subdivided  as  follows! 

x j = a ■»  jh  , j^O.l , . . . ,N 

h • (b-a)/(N+l) 

We  call  the  points  the  mesh • 


An  approximation  to  the  solution  of  (4. 1.1), (4. 1.2) 
can  now  be  obtained  by  replacing  these  expressions  with  a 
set  of  finite-difference  equations.  The  first  such  method, 
Stormer’s  method  (see  Henricil_23j#  Chapter  7),  results  in 
finite-difference  equations  of  the  form 

- Yi_i  +2Yi  " Yi  + l 4 h2(q(Xi)  - A)Y.  = 0 (4.1.3) 

i=l, 2, . . . ,N 

Y0  = Vl  = °* 

The  discretized  function,  {.Y^},  i6  an  approximation 
to  some  eigenfunction  of  (4.1.1),  and  A,  an  approximation 
to  the  corresponding  eigenvalue,  (see  KellerL25].  Chapter 
5,  for  a slightly  more  general  discussion) 

If  we  define  Y to  be  the  N-dimensional  column  vector 


Y = (Y^,  Yg.  ....  Yjj ) 


and  A to  be  the  N x N matrix  A = (a. .)  where 

^ J 


aj.M  = aj.M  = _1 


e.j  = 2 + h q(x j) 
a^  = 0 if  |i-j|  > 1 

then  (4.1.3)  can  be  written  as  the  following  matrix  problem. 


AY  - hzAY  = 0 


(4.1.4) 


The  matrix  A has  the  useful  properties  of  being  both  sym- 
metric and  tridiagonal. 


4.1.1  Truncation  Error 
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The  local  truncation  error  of  the  method  (4.1.3)  at 

x ^ is  defined  as  the  amount  by  which  the  exact  solution, 

(y,M , fails  to  satisfy  the  discretization,  i.e. 

T(xi)  = -y(xi_1)  + 2y(xi)  - y(xi+1)  ♦ h2tq(x)  - X]y(xi) 

= -y(xi_1)  ♦ 2y(xi)  - y(xi+1)  + h2y  (x.^ 

(4.1.5) 

Expanding  yfx^^)  and  y(xi+1)  in  powers  of  h we  find 

that 

“C(Xi)  = - J2  hViV)(0)’  xi_1  < 0 < xiU 

i=l , 2,  • • • 

If  we  now  define  Tj^y)  = ~2  (xi ) - yg  h2y^iv^(0) 

and  t(y)  = (T^y) , X 2(y) . • • • * tN(y)  we  can  write  (4.1.5) 

as  the  following  matrix  equation. 

Ay  - h2Xy  = h2  f (y)  (4.1.6) 

y(x0)  - y(xNfl)  - 0 

where  y = (y(x1),  y(x2) y(xN) 

4.1.2  Convergence 

The  matrix  eigenvalue  problem,  (4.1.4),  possesses 
exactly  N eigenvalues  while  the  differential  equation, 
(4.1.1),  has  an  infinite  number.  The  problem  then  is  de- 
termining which  of  the  exact  (differential  equation)  eigen- 
values are  being  approximated,  and  how  accurate  that  approx- 
imation is. 
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Keller  provides  us  with  a theorem  which  clarifies 
this  point.  We  present  it  here  without  proof.  The  vector 
norm  used  is  the  2-norm. 


THEOREM  4 . 1 » (Keller)  For  each  fixed  eigenvalue  X 


of  (4.1.1)  with  corresponding  eigenfunction  y(x),  there 


exists  an  eigenvalue  h v\  of  A such  that 

I/S. 


K | < 

" llyll 


It  is  an  easy  matter  now  (see  for  example  Keller 
t>.  134-13e0  to  show  that 


IA  - X|  < 0(h4) 


Thus  we  see  that  any  exact  eigenvalue  can  be  approximated 
to  a desired  accuracy  by  taking  h small  enough. 


4.1.3  Solution  of  the  Matrix  Problem 

If  one  wants  to  compute  all  the  eigenvalues  and 
eigenvectors  of  a matrix,  QR  methods  have  become  the  pre- 
ferred technique.  Stewart[.39j  gives  a thorough  treatment 
of  both  the  explicitly  shifted  and  the  implicitly  shifted 
algorithms . 

It  is  shown,  however,  in  Chapter  that  during  the 
adaptive  solution  process  the  final  matrix  problem  to  be 
solved  depends  on  the  particular  eigenvalue  under  consid- 
eration. In  general,  the  final  mesh  and  thus  the  final 
matrix  problem  will  be  different  for  each  eigenvalue. 

An  algorithm  such  as  QR  is  therefore  mucSi  too 
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expensive  to  be  used  for  our  purposes  since  we  solve  for 
only  one  eigenvalue  from  any  particular  matrix  problem. 

Instead  we  use  a technique  called  Rayleigh  quotient 
iteration.  An  exhaustive  analysis  of  this  technique  appears 
in  a series  of  papers  by  0strowskii.3l]  • We  refer  to 
results  from  these  papers  frequently  throughout  this 
chapter. 

The  algorithm  we  use  is  as  follows. 

ALGORITHM  4.2>  Let  A be  a real,  symmetric  matrix 
and  let  AQ  be  an  approximation  to  A.  (Discussion  of  how 
to  set  such  an  approximation  is  deferred  until  Chapter  6.) 
Let  € be  the  desired  tolerance  level. 

1.  Let  Yq  = (1,  1 1)*  [.arbitrary] 

2.  i = 0 

3.  Solve  (A  - h2AiI)  ui+1  = Yi  for  ui+1 

4.  Compute  Ai+1  = uj+1  Au.^  / h2u f +1ui+1 
(Rayleigh  quotient) 

5.  Compute  Y.+1  = u.^  / (u*+1  ui+1)* 

(Normalization  of  the  eigenvector) 

6.  If  IAi+1  -AJ  < * STOP  else  i = i+  1 and 
GO  To  3. 

Ostrowski  shows  that  if  A is  symmetric  this  process 
converges  cubically  to  an  eigenvalue,  i.e. 

IA1+1  -A|  . clAi  - Al3 

The  shooting  methods  described  in  Chapter  3 rely  on 
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Newton's  method  for  convergence  and  thus  exhibit  only 
quadratic  convergence  to  the  eigenvalue. 

For  the  remainder  of  this  chapter  we  assume  that  the 
approximate  eigenvalue  required  to  start  Algorithm  4.2  and 
all  other  algorithms  to  follow  is  always  available.  In 
Chapter  6 we  present  an  algorithm  for  efficiently  computing 
a good  approximation. 

4.1.4  Error  Estimate  and  Deferred  Correction 

theorem  4.1  showed  that  by  taking  h small  enough  we 
can  get  A arbitrarily  close  to  a.  Additionally,  for  any 
particular  h we  are  able  to  estimate  the  error  in  A. 

Recall ing  (4.1.6), 

Ay  - h2Ay  = h2  t(v)  (4.1.6) 

we  let  v - Y + -'y  and  a A + . Substituting  in^o 

(4.1.6)  and  multiplying  by  Y*,  we  obtain 

YflAY  - h2AY»  + LA^Y  - h^Yl^y 

- h^AY^y  * h2Yt  X(v)  (4.1.7) 

The  first  term  of  (4.1.7)  is  equal  to  zero  by 
(4.1.4).  "’he  second  term  is  also  zero  since  A is  sym- 


metric . 
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Solving  for  we  get 

- - Y*  t(v)  / Y^y 

which  is  an  exact  expression  for  the  error  in  A.  This 
expression  is  not  a practical  one,  however,  because  it  in- 
volves v,  the  exact  solution.  If  in  (4.1.7)  we  assume  that 
^Aay  is  negligible  then  an  approximation  for  becomes 

* - yt  t(y)  / Y*Y  (4.1.R) 

The  vector  t(y)  is,  of  course,  also  a function  of 
the  exact  solution  y,  but  TT(y)  , too,  can  be  approximated 
by  taking  the  appropriate  differences  of  the  C Y ^ } or  equiv- 

M 

alently  the  lY^|.  A method  for  approximating  the  trunca- 
tion error  is  presented  in  Chapter  5*  For  now  we  assume 
only  that  we  have  available  a metnod  which  produces  the 
approximation 

Vy)  = T i ( y ) £ 1 ♦ 0(hr)  ' , r --  1. 

Let  us  define  T(Y ) - (T^Y).  T2(Y) T^( Y)  )t. 

We  can  now  introduce  an  algorithm  known  as  deferred  cor- 
rection and  made  popular  by  FoxLl8].  It  enables  us  to 
improve  upon  the  matrix  eigenvalue  and  eigenvector  and  is 
described  as  follows. 
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ALGORITHM  4,3«  Let  Y and  J\  be  the  solution  to  the 
matrix  problem  (4.1.4). 

1.  Estimate  T(y)  by  computing  T(Y). 

2.  Estimate  aA  by  computing  aX-  - Y*T(Y)  / Y*Y. 

3.  Let  A ^ A ♦ aT. 

4.  Solve  (A  - h2Vl)  y « T(Y)  for  y. 

The  solution  (y.T)  is  now  presumed  to  be  a better 
approximation  to  (y,A).  The  work  by  FischerLl6],  although 
it  is  concerned  specifically  with  a higher  order  method, 
implies  that  if  Y were  found  to  be  equal  to  y * 0(h  ),  then 
y = y + O(h^).  From  what  follows  Theorem  4.1  we  know  aA 
to  be  0(h).  Thus,  one  correction  by  Algorithm  4.3  also 
yields  V = A ♦ 0(h**)  . 

Algorithm  4.3  can  also  be  used  iteratively  by  re- 
computing aA  and  T(y)  and  repeating*  however,  Fox  argues 
that  usually  only  one  correction  is  necessary  to  get  high- 
ly accurate  results. 

Numerical  results  using  the  Stormer  method  are  pre- 
sented at  the  end  of  Section  4.2,  at  which  time  they  are 
compared  with  the  higher  order  Numerov  method  of  that 
section. 

4.2  The  Numerov  Uniform  Mesh  Method 

As  the  next  step  toward  our  adaptive  method  we 
introduce  a higher  order  approximation  to  (4.1.1),  (4.1.2). 
While  this  method  still  uses  a uniform  mesh  it  possesses 
some  properties  different  from  those  of  the  Stormer  method 
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and  thus  serves  to  introduce  some  concepts  essential  to  the 
adaptive  method. 

The  discretization  we  use  is  as  follows. 

- Yi-1  * 2Yi  - Ym  * * l0Vi  * 


fi+lYi-a)/12  = °*  i=l , 2, . . . ,N 


Yo  a YN+1  = 0 


(4.2.1) 


where  = (q^  - A).  This  is  known  as  the  Numerov  method, 
When  we  write  (4.2.1)  in  matrix  form  we  have 


AY  - A BY  " 0 


where  Y is  as  before,  A is  the  N x N matrix,  A = 

a j . 3-1  = -1  + h2^i-l  / 12 

aj,j+l  = _1  + h2qi+l  / 12 
a j j * 2 10h2qi  / 12 

ai<5  = 0 if  |i-j|  > 1 

and  B is  the  N x N matrix,  B ^ ^ ) , with 

■ bJ.M  ’ h2/12 

bjj  = 10h2/l2 

b.j  =0  if  | i - j | > 1. 


(4.2.2) 


(a^),  with 


The  matrix  B is  tridiagonal,  symmetric,  and  positive- 
definite,  but,  although  A is  still  tridiagonal,  we  see  that 
it  is  no  longer  symmetric.  (Figure  4.1) 
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A = 


1 

12 


24tlOh  qx 
-12+h2q^ 


24+10h2q2  -12+h2q^ 

-12+h2q2  24-»-10h2q3  - 


^4 


Figure  4.1 
The  Matrix  A 


Matrix  A possesses,  however,  another  quite  useful 
property.  First,  we  make  the  following  definition. 


DEFINITION » (Wilkineon[42])  Let  C = (c^)  be  an 
n x n,  real, tridiagonal  matrix  defined  by  = a^, 

ci,i+l  = ^i+1*  ci*l, i = Yi+1’  cij  = 0 otherw*se. 

Then  C is  auasi-symmetric  if  t* > 0 for  i=l,2,..,n. 


A quasi -symmetric  matrix  can  be  transformed  into  a 
real,  symmetric  matrix  as  follows. 

Let  dj^  = 1*  t YiY2*  * * Yi^  i®  2*  * i ^ij  = ^ 

for  i / j.  Then  D - <d14)  is  such  that  D”^CD  = T is  sym- 
metric  with  tu  = «j,  t1(ln  = tJ<ia  = (6lnYln>}. 

We  can  now  make  the  following  observation  about  the 
matrix  A. 


THEOREM  4,4 i For  h sufficiently  small,  A is  quasi- 
symmetric. 

Proof i For  the  theorem  to  be  true  we  must  find  an 
hQ  such  that  for  all  h < hQ  we  have 


i 1 f 


r 


*ravm 


i 
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ai,i+l*ai+l,i  3 (_1  * *2qin/l2)(-l  + h2q./l2)  > 0. 

(4.2.3) 

If  we  let  Q = max  |q.|  then  by  choosing  h_  such  that 

J J ° 

-1  ♦ h2Q/12  <0  or  hQ  < (12/Q)^ 

each  term  on  the  right-hand  side  of  (4.2.3)  will  be  neg- 
ative for  h < hQ  and  the  theorem  is  proved. 

One  further  property  of  the  matrix  problem  (4.2.2) 
which  we  must  show  is 

PROPERTY  I t For  h sufficiently  small  the  matrix 
problem  (4.2.2)  has  N real,  distinct  eigenvalues. 

Proof » The  reader  is  referred  to  Corollary  4.18 
of  Theorem  4.16  for  the  proof. 

Because  the  Numerov  method  is  merely  a special  case 
of  the  non-uniform  mesh  method  presented  in  Section  4.3, 
several  results  thoughout  the  remainder  of  this  section  are 
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Defining  t(y)  = (t (x^) ,1(x2) , . . . ,T(xN)  )t  we  can 

write 

Ay  - A By  - f(y)  (4.2.4) 

4.2.2  Convergence 

The  task  of  showing  that  the  matrix  eigenvalues 
converge  to  the  exact  eigenvalues  is  slightly  complicated 
by  the  quasi -symmetry  of  A.  Kellerl.25]  presents  a theorem 
which  states  that* 

If  A is  symmetric  and  B symmetric-positive-definite 
then  for  each  eigenvalue  X of  (4.1.1)  and  corresponding 
normalized  eigenfunction  y(x)  there  exists  an  eigen- 
value A of  B-1A  such  that 

I A - a | < ||B“15J*ll  t(y)ll  / llyll 

We  have  shown,  however,  that  the  Numerov  method 
produces  an  A which  is  quasi -symme trie . The  problem  of 
adapting  Keller's  theorem  to  accomodate  the  quasi -symmetry 
of  A appears  quite  difficult.  Instead  we  prove  convergence 
in  a slightly  different  manner. 

The  statement  of  our  result  is  as  follows* 

PROPERTY  lit  For  the  Numerov  method,  (4.2.1), 
there  exists  , for  every  eigenvalue  X and  corresponding 
eigenfunction  y(x)  of  the  differential  equation,  an  eig- 
envalue A of  the  matrix  problem  (4.2.2)  and  a correspond- 
ing eigenvector  Y such  that 

(A  - X)  = 0(h4) 

and  Y^  = y(x^)  ♦ 0(h4),  i=0,l, . . . ,N+1. 
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Proof » See  Corollary  4.20. 

We  see  now,  as  with  the  Stormer  method,  that  by 
taking  h sufficiently  small  we  can  approximate  the  exact 
solution  to  within  an  arbitrary  tolerance. 


4.2.3  Solution  of  the  Matrix  Problem 

As  with  the  Stormer  method  of  Section  4.1,  we  want 
to  solve  the  matrix  problem  (4.2.2)  for  only  one  eigenvalue 
and  eigenvector  at  a time.  We  again  use  Rayleigh  quotient 
iteration. 

It  is  important  to  note  here  that  the  fact  that  our 
problem  has  real  eigenvalues  is  essential  to  us  now. 
Rayleigh  quotient  iteration,  as  defined  here,  can  only 
converge  to  a real  eigenvalue.  Although  it  was  not  men- 
tioned at  the  time,  the  real  eigenvalue  property  also  holds 
for  the  Stormer  method  s5nce  all  eigenvalues  of  a sym- 
metric matrix  are  real. 

We  present  now  the  algorithm  for  solving  the  gener- 
alized matrix  eigenvalue  problem  (4.2.2)  by  Rayleigh  quo- 
tient iteration. 


ALGORITHM  4.S:  Let  A and  B be  the  real,  N x N 
matrices  described  above,  Let  A q be  an  approximation  to 

1.  Let  BYq  = (1,1, ... ,1) 1 [arbitrary] 

2.  i = 0 


3.  Solve  (A  - A^lu^j  - BY^  for  u 

4.  Compute  ~ ^i4i^^*i+i  ^ ^i+1 


i*l 

Bu 


i*l 
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5.  Compute  ^i+l  = ^i+1  ^ 

6.  If  lAj+i  “ Ail  < e STOP  else  i = i+1  and 
GO  TO  3. 

Ostrowski  also  discusses  this  process  in  his  series 
and  shows  that, for  A symmetric  and  B symmetric-positive- 
definite,  cubic  convergence  of  J\  ^ to  A is  again  present. 
Unfortunately  this  result  does  not  immediately  carry  over 
to  the  case  in  which  A is  quasi -symmetric . However,  we 
show  in  Section  4.3.6  that  because  of  the  interdependence 
between  the  matrix  problem  and  the  differential  equation, 
it  is  still  possible  to  obtain  cubic  convergence  under 
certain  conditions. 

4.2.4  Error  Estimate  and  Deferred  Correction 

The  derivation  of  an  error  expression  for  (A  - ) 

is  again  possible. 

Starting  with  (4.2.4), 

Ay  - KBy  * ?(y),  (4.2.4) 

letting  y = Y Ay  and  i » A 4 ^ . and  substituting  into 
(4.2.4)  while  multiplying  by  Y*  we  get 

Y*A(Y  ♦ Ay)  - (A  + AX)YtB(Y  ♦ Ay)  - Y*  t(y)  . 

Multiplying  out,  we  have 

Y^AY  - ABY)  ♦ Ytt  AAy  - /\Bby)  - 

- AXYtBY  - AiY^BAy  - Yt  T(y)  (4.2.5) 
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The  first  term  is  zero  directly  from  (4.2.2).  The 
fourth  term  contains  the  product  and  so  can  be  consid- 

ered negligible  in  relation  to  the  other  terms. 

The  second  term  we  can  also  consider  to  be  ne^l i - 
gible  by  the  following  reasoning. 


AY  = A BY 

YtAt  - AYtBt  = A YtB 

YtAtAy  - AyVv 

YtMy  - AyVv  - Yt/Uy  - YtAt^y 
Yt(A  - AB)^y  Yt(A  - At),'^y 

However,  the  matrix  (A  - A^)  consists  of  0*s  on  the  diag- 
onal  and  terms  like  h (q^  - q^^/l?  on  the  super  and  sub- 
diagonals . 

With  q(x)  being  continuous  on  these  terms 

are  0(h-^)  . Thus,  when  we  solve  (4.2.*)  for  -^a , eliminating 
the  zero  and  negligible  terms,  we  have 

« . **  iiil  . ?.*.(*  -..  **)■** 

YtBY  YfBY 

and  we  find  that  the  first  term  is  Oth^)  while  the  second 
is  0(h<’).  This  leads  us  to  the  following  error  estimate. 

•AA  * - Y1  t(y)  / YtBY . (4.2.6) 

We  note  again  that  to  use  this  estimate  we  must 
approximate  t(y) , and  such  an  approximation  is  presented 


■■■ 
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in  Chapter  5«  For  now  we  assume  that  T^(Y)  =T^(y)* 

(l  ♦ 0 (hr) } • r > 1,  is  available. 

The  error  estimate  (4.2.6)  now  makes  available  the 
means  for  correcting  the  Numerov  method  in  the  fashion  of 
Algorithm  4.3. 

ALGORITHM  4.6 » Let  Y and  J\  be  the  solution  to  the 
matrix  problem  (4.2.2). 

1.  Estimate  T(y)  by  computing  T(Y) . 

2.  Estimate  aX.  by  computing  aT  - -Y^TfY)  / Y^BY . 

3.  Let  T = A * a*  . 

4.  Solve  (A  - XB)y  - T(Y)  for  y. 

The  improvement  to  the  eigenvector  by  making  one 
such  correction  was  shown  by  FischerLl6]  to  be  from  an 
error  of  O(h^)  to  O(h^) . A similar  improvement  is  ex- 
perienced in  the  eigenvalue.  (See  results  below) 

4.2.5  Numerical  Results 

In  Tables  4.1  and  4.2  are  presented  some  numerical 
results  for  the  following  problem,  known  as  Weber's 
equation. 

y ♦ (1  - x2)y  = 0 
y(0)  = y(l)  = 0 

Table  4.1  shows  results  for  XQ  « 10.1511640305 
while  Table  4.2  represents  * 89.1543424562. 

Each  table  shows  the  following  for  N = 8,16,32, 

64,  and  128. 


r 


4C 


F 

1.  Absolute  value  of  error  estimate  (4.1.8) 

2.  I A * M for  Stormer's  method  (E  ) 

s 

3.  ESA2 

4.  Absolute  value  of  (4.2.6) 

5.  I A - A I for  the  Numerov  method  (E  ) 

n 

6.  En/h4 

7.  |X  - X|  for  Stormer's  method  (E„  ) 

sc 

(Note i this  is  the  corrected  eigenvalue.) 

8.  Et/ 

9.  1^  - 4 | for  Numerov  method  (E_) 

nc 

1°.  Enc/h6 


Table  4.1 


Results  for 


N 

14.1.81 

Es 

Es/h2 

Esc 

E /h4 
sc 

8 

1.237D-1 

1.263D-1 

8.083 

2.654D-3 

10.87 

16 

3.152D-2 

3.169D-2 

8.113 

1.675D-4 

10.98 

32 

7.920D-3 

7.929D-3 

8.119 

1.049D-5 

11.00 

64 

1.982D-3 

1.983D-d 

8.122 

6.56ID-7 

11.01 

128 

4.957D-4 

4.957D-4 

8.123 

4.11  D-8 

11.03 

N 

14.2.61 

En 

Enc 

En</h< 

8 

1.122D-3 

1.066D-3 

4.366 

5.541D-5 

14.53 

16 

6.596D-5 

6.628D-5 

4.344 

3.I66D-7 

5.312 

32 

4.130D-6 

4.140D-6 

4.341 

6.3 

D-9 

6.765 

64 

2.5840-7 

2.600D-7 

4.362 

1.0 

D-10 

6.872 

128 

1.615D-8 

2.000D-8 

5.369 

5.0 

D-ll 

*2199. 

* xo 

correct  to 

all  available  digits 

i 


Table  4.2 
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Results  for 


N 

14.1.8! 

£ 

s 

Vh 

E „ 
sc 

4 

\ 

O 

to 

W 

8 

8.130DO 

9.811D0 

627.9 

1.681D0 

6885. 

16 

2.434D0 

2.539D0 

650.0 

1.681D-1 

7556. 

32 

6.329D-1 

6.403D-1 

655.7 

7. 37 2D- 3 

7731. 

64 

1.599D-1 

1.604D-1 

657.0 

4.634D-4 

7774. 

128 

4.010D-2 

4. 01 2D- 2 

657.3 

2.901D-5 

7787. 

N 

14.2.61 

E 

n 

E 

nc 

En</h6 

8 

6.992D-1 

7.488D-1 

3067. 

4.962D-2 

13008. 

16 

4.877D-2 

4.521D-2 

2963. 

3.562D-3 

59760. 

32 

2.785D-3 

2.797D-3 

2920. 

1.264D-5 

13572. 

64 

1.740D-4 

1.744D-3 

2926. 

4.299D-7 

29543. 

128 

1.089D-5 

1.089D-5 

2923. 

7.2  D-9 

31666. 

Tables  4.1  and  4.2  illustrate  several  points.  We 

2 4 

see  the  improvement  in  order  from  0(h)  to  0(h  ) when  the 
Numerov  method  is  used.  We  note  that  both  methods  experi- 
ence the  expected  increase  in  order  when  the  correction  is 
made.  Finally,  we  note  that  the  error  estimates  do  appear 
to  approach  asymptotically  the  actual  error. 

4.3  A Non-Uniform  Finite-Difference  Method 

It  is  the  contention  of  this  thesis  that  a more 
efficient  method  of  solving  the  S-L  problem  is  to  use 
finite-difference  methods  adaptively  with  a non-uniform 
mesh.  In  preparation  for  such  a method  we  present  a non- 
uniform  discretization  and  describe  its  application  to 
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our  problem.  Our  selection  of  the  particular  form  of  the 
discretization  was  influenced  by  the  efficiency  of  the 
Nuraerov  method  of  Section  4.2,  We  wish  to  retain  the  tri- 
diagonal property  of  the  matrices,  the  cubic  convergence 
of  Rayleigh  quotient  iteration,  and  the  O(h^)  behavior 
of  the  error.  As  we  show  in  this  section,  the  following 
form  does  preserve  these  properties. 

Let  n = be  a mesh  defined  on  La,b]  with 

a = x0  < -c  ...  < xN  < xN41  b.  Let  h.  x.  - x.^, 
i=l,  2 N+l . 

We  now  introduce  the  following  discretization  of 
the  differential  equation  (4.1.1). 


“01*1-1  * 2*i  * “21*1*1  * d0i*  l-l  * “’ll*  1 


* S2iI  ln  . 0.  1-1.2 N (*.3.1> 


where  a2^,  3Qi,  3^,  and  32i  are  real  coefficients 

•• 

to  be  determined  and  - (q(x^)  - A)Y^. 

If  we  now  require  that  this  approximation  be  exact 
for  y(x)  - 1,  x,  x , xJ , and  x we  get  the  following  set 
of  equations « 


y(x) 

y(x) 

y(x) 


1,  y (x)  * o 

“01  * 2 * “21 
#• 

X,  y (x)  =■  0 

“01*1-1  * 2xl 
*2.  y'(x)  * 2 


0 


a2iXi+l  ' 0 

I 

“21xl*l  *»01  ‘ “li  * 29 21 


0 


y(x)  = x3,  y (x)  = 6x 
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aOixi-l  * 2xi  * a2lxi*l  + 680ixi-l  * 

<692lxi*l  ' 0 
y(x)  - x1*,  y (x)  = 12x2 

aOixi-l  * 2xi  + a2lxi»l  * 128oixi-l  - 128lixi 
* l»2ixiU  * 0 

Solving  the  five  equations  simultaneously  and  re- 
placing the  x^'s  with  appropriate  h^*s,  we  get 


a0i  = "2hi+l  / ^hi  * hi+i) 

a2i  = ’2hi  / (hi  * »W 

0Oi  = ~hi+l(hi-*-l  " hihiU  ' hi}  ^ 6(hi  4 hi«l* 

Bli  = (hi  + 3hihi+l  ♦ h2+1)  / 6 

0 2i  = 'hi(hi  ‘ hihi  + l " hit-l)  / 6(hi  4 hi-fl) 


Thus  the  difference  equation  is,  in  general,  different  at 
each  point  x^. 

These  coefficients  are  similar  to  those  used  by 
BrownL?]  to  solve  non-linear  boundary  value  problems.  He 
does  not,  however,  consider  eigenvalue  problems. 

As  with  the  previous  discretizations  we  can  now 
write  equations  (4.3,1)  in  matrix  form  to  get 

- 

(A  - AB)  Y » 0 (4.3.2) 

where  A = (a.  .)  and  B -•  (b.  .)  are  N x N real,  tridiagonal 
* J ^ J 

matrices  and  Y = (Y^)  is  a real  column  vector  of  length  N 
with 


II 

D 


The  problem  now  becomes  that  of  finding  the  eigenvalues, 

, and  eigenvectors,  Y,  of  the  matrix  problem  (4.3.2). 

We  note  that  the  symmetry  and  positive -definite 
properties  of  A and  B are  no  longer  present.  Only  the  tri- 
diagonal property  remains.  If  some  further  simplifying 
properties  were  not  available  for  A and  B,  a general  algo- 
rithm such  as  the  QZ  algorithm  (see  StewartC39],  Chapter  ?) 
would  be  required  to  solve  the  problem.  We  are,  however, 
able  to  show  several  interesting  properties  of  the  matrix 
problem,  each  of  which  serves  to  simplify  the  solution 
process. 

4.3.I  Properties  of  the  Matrices  A and  B 

First,  we  note  two  obvious  facts  about  A and  B. 

Since  h^  is  always  positive  it  follows  that« 


PROPERTY  4.7»  and  a2^  are  negative  for  all  i. 
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PROPERTY  4.8 1 0^  is  positive  for  all  i. 

These  properties  are  immediately  obtained  from  the  expres- 
sions for  a 0^,  a2^,  and  *li- 

In  the  paper  by  Brown,  mentioned  above,  he  requires 
that  and  *2i  always  be  positive.  This  is  essential 
to  his  proof  of  the  convergence  of  Newton's  method  for 
solving  his  non-linear  problems. 

With  the  following  lemma  we  show  the  severe  restric- 
tions this  places  on  the  relative  sizes  of  adjacent  inter- 
vals, h^  and  h.^. 

Let  h.^  = a*h^  where  a > 0. 

LEMMA  4.9» 

i)  dQi  - 8 2i  ~ 1/12  when  h^  = h^+1  (a  = 1). 
ii)  Both  and  d ^ are  > 0 if  and  only  if 

V5  ~ 1 < a < V3  1 1 
2 a 2 

m)  0oi02i  > 0 if  and  only  if  dQ^  > 0 and  d2j  > 0. 

Proof i 

i)  This  is  immediately  obvious  by  substituting 
hi  = hi+l  *nt0  the  expressions  for  dQ^  and  d2^.  Also 
here  we  note  that  0^  * 10/12  when  tv  = h^4^.  These  coef- 
ficients are  of  course  what  we  expect  since  when  h^  • h.^ 
the  discretization  is  simply  the  Numerov  method  of  Section 
4.2. 


ii)  d, 


' -hin<hi.1  - Vin  - hi>  / 6<hi  • hi-i> 

(-a2  ♦ a * l)h.n  / 6(h.  ♦ h-M)h? 


So  > 0 if  and  only  if  -a  ♦ a * 1 > 0 which 
occurs  if  and  only  if 


< a < 


, V5_L_1 


2 


Similarly  ^2i  > 0 2111(1  only  i** 


a < ~^-S. ~ .X  or 


Thus,  both  and  t)2.  > 0 if  and  only  if 


< a < 


V5_|_l 


.618  < a < 1.618  (approximately). 


iii)  By  steps  similar  to  those  of  ii)  we  find  that 


< 0 if  and  only  if 


_L_i  or 

2 °r  2 


and  that  < 0 if  and  only  if 


. 

There  is  no  positive  value  of  a which  can  satisfy 
both  and  so  we  conclude  that  and  o ^ cannot  both  be 
negative . 

iii)  follows  from  ii). 


The  significance  of  this  result  is  as  follows.  If 
we  were  to  require,  as  Brown  did,  that  the  entries  of  B all 
be  positive,  then  we  must  restrict  the  ratio  of  adjacent 
steps  to  that  of  ii).  This  is  so  restrictive  that  steps 
cannot  even  differ  by  a factor  of  2. 

If  we  lift  the  restriction  and  require  only  that 
steps  be  changed  by  an  integer  factor  (i.e.  a = n or  1/n, 
n a positive  integer)  we  get  the  following  corollary. 

COROLLARY  4.10»  If  adjacent  step  sizes  always  dif- 
fer by  an  integer  factor  then  either  2 ^ - 1/12  or 

B0i82i  < 

We  see  in  Section  5.3  that  this  condition  is  auto- 
matically met  with  our  method  of  refining  the  mesh. 


4. 3. 1.1  Quasi -symmetry  of  A 


With  the  above  properties  available  it  is  now  pos- 
sible to  show  that  under  certain  conditions  the  matrix  A 
is  quasi -symmetric,  as  was  the  case  with  the  Numerov  method. 

Letting  h = max  h.  we  have  the  following 1 
l<i<N+l  1 


: 

1 


THEOREM  4.11 i For  h sufficiently  small,  A is  quasi- 
symmetric. 

Proof » For  A to  be  quasi -symmetric  we  need  to  find 
an  hQ  such  that  for  all  h < hQ 

2.1-1  *e2.1-l,'l)  * 0 for  a11  l' 


(a01  * B01<'i-l)(,, 
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Let  h^  = ph  and  h.^  * rh  where  0 < p < 1 and 
0 < r < 1. 

We  first  show  that  for  h sufficiently  small 

“oi  * "ol'l-l  < 0 or  -“oi- 

From  Property  4,7  it  is  known  that  < 0,  We 
also  note  that 

30i  = "hi+i^hi+i  * hihi+l  “ hi^  / + hi+i^ 

‘ “oi(hi.i  - hihi-i  - hi>  / 12 
- -aQi(h2/l2) (p2  ♦ pr  - r2) 

So  we  want 

-aQi(h2/l2) (p2  + pr  - r2)qi_1  < -aQi 
or  dividing  by  -aQ^/l2t  which  is  positive,  we  get 

h2(p2  ♦ pr  - r2)qi_1  < 12.  (4.3.3) 

Letting  z(p,r)  - p * pr  - r , we  wish  to  find 

max  |z(p,r)|  where  R * [(p.r)  I 0 < p < 1,  0 < r < l] . 

R 

The  critical  point  occurs  where 
dz/dp  = dz/«tr  * 0 


or  where 


2p  + r s 0 
p - 2r  * 0 

or  at  (p,r)  • (0,0).  (z  « 0) 

Hence,  the  maximum  occurs  on  a boundary,  namely  at 


r 
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(p,r)  - (l.i)  where  z ~ 5/4. 

max  | z ( p , r ) | - 5/4 . 

R 

If  we  also  define  Q - max  |q(x)|,  we  find  that 

a<x<b 

(4.3.3)  holds  if 

h < (48/5U)*. 

We  also  need  02  j^q^  < -a2i  whi-ch  leads  to  the  same 
restriction  on  h. 

Thus  letting  hQ  = (48/5Q)^  we  see  that  for  h < hQ 
A is  quasi - symme tr ic . 

4. 3. 1.2  Properties  of  B 

We  have  already  noted  that  the  matrix  B is  no  longer 
symmetric,  positive-definite,  however,  B does  still  possess 
several  of  the  important  properties  of  a positive-definite 
matrix  and  we  present  them  here. 

The  first  such  property,  although  important  in  it- 
self, is  mainly  used  to  prove  several  facts,  and  so  is 
introduced  here  as  a lemma. 

LEMMA  4. 12 i B is  strictly  diagonally  dominant. 

Proof t We  must  show  that 

|t*li  I » l0Qil  + l02il  for  i-1,2 N. 

Ii*0i1  4 l02i'  * ^hi+l  * hihi+l  " hihi+l*  + 

|hi  “ hihi*l  " hihin,J/6(hi  * hi*l) 
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< [h?  + 2hihi+1  + 2hihi+l  4 hi+l  J/^hi  4 hi+l^ 

< [h^  + 3h?hi+1  + 3hj.h2+1  + hi+iV6(h^  ♦ h^j) 

= L(b^  ♦ h^.^)^]  / 6(h^  + h^+^) 

= [h?  ♦ -2hihi+1  + h2+1]  / 6 

< [hf  ♦ 3hihi+1  + h?+1]  / 6 

* ~ t ^li ^ * 

The  next  result  we  need  for  our  ultimate  statement 
about  the  eigenvalues  of  (4.3*2)  is  one  about  the  de- 
terminants of  B and  its  leading  principal  sub-matrices. 

The  desired  result  appears  as  a corollary  to  the  following 
theorem. 

THEOREM  4.13i  If  B is  a real  n x n matrix  and  is 
strictly  diagonally  dominant  with  positive  diagonal  ele- 
ments. then  det(B)  > 0. 

Proof i By  a theorem  of  Gerschgorin,  taken  here 
from  Varga{.4l],  we  know  that  the  eigenvalues,  [X^}#  of  a 
strictly  diagonally  dominant  matrix  with  positive  diagonal 
entries  satisfy 

Re(k^)  > 0,  1 < i < n. 

This,  of  course,  implies  that  all  real  eigenvalues  are 
positive . 


Also  it  is  known  (see  for  example  StewartC  39  ] , 
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p.  269)  that  the  complex  eigenvalues  of  a real  matrix  occur 
in  conjugate  pairs  and  that  the  product  of  a number  and  its 
complex  conjugate  is  a positive  real  number. 

Thus, 

" X , > 0.  (4.3.4) 

j=l  3 

But 

n 

it  X . - det(B) 

j=l  3 

as  the  following  argument  shows. 

Let  p(X)  be  the  characteristic  polynomial  of  B, 
i.e.  p(X)  = det(B  - XI).  We  know  that  p(X)  can  be  factored 
into 

p(X ) *-=  ( -l)n(X  - X1)(X  - X2)  ...  (X  - xn) 

« det  (B  - XI) 

Letting  X = 0 in  this  identity,  we  get 

<-l)n(-l)n  I X = det  (B).  (4.3.5) 

j*l  3 

And  so  combining  (4.3.4)  and  (4.3.5)  we  see  that 
det  (B)  > 0. 

This  theorem  can  also  be  shown  to  be  an  immediate 
corollary  to  the  following  theorem  by  G . B.  Priced 34], 
stated  here  without  proof. 


THEOREM  4.14 i Let  A = 


real  elements  such  that 
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(a^.)  be  an  n x n matrix  with 


then 


ii 


n 

> £ | a 


J-l 

j/i 


U1 


0 ^ • • • ni^  < det  (A)  < MXM2 . . .Mn 


where 


n 


m.  - | a.. I - £ |a.. 

1 11  j-i+1  1J 

n 

M.  --  I a.  . | ♦ £ la.. 


ii 


j-in  13 

If  we  now  define  to  be  the  leading  principal 
submatrix  of  B of  order  i,  we  get  the  following: 


COROLLARY  4.15:  For  the  matrix  problem  (4.3.2) 
det  ( ) > 0,  i=l,2,...,N. 

Proof:  Each  B^  is  a real  i x i matrix.  From  Lemma 
4.12,  each  Bi  is  strictly  diagonally  dominant.  From 
Property  4.8,  each  B^  has  positive  diagonal  elements. 

By  Theorem  4.13,  det  (B^)  > 0,  i-1,2 N. 

It  is  possible  to  develop  a 3-term  recurrence  re- 
lationship for  det  ( B^ ) . Because  such  a relationship  is 
needed  in  a later  proof,  we  introduce  it  here. 

Let,  us  define  det(BQ)  - 1 and  det(B^)  --  8^. 

Recalling  that  B,  and  thus  B^,  are  tridiagonal  we 
see  that  det(B^)  can  be  computed  by  expanding  B^  about  row 


i to  yield 
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det( B^)  ■ 9lidet(Bi_1)  - 0QiP2  j^jdetlBj,  _2) 

i-2,3 N (4.3*6) 

Expression  (4.3.6)  is  the  3-term  recurrence  relationship 
we  need. 

4.3.2  Truncation  Error 

The  next  property  of  our  discretization  we  wish  to 
look  at  is  the  truncation  error.  It  can  be  computed  by 
first  noting  that 

T(xi)  - «oiy(xi_1)  ♦ 2y ( x i ) ♦ a2iy(xi-^l)  + 


h *i  •• 


*0iy  (xi-l 

) ♦ P 

liy  (x 

i>  ♦ *21* 

(xi.i' 

(4.3.7) 

-2hi+1y(x. 

‘ hi 

)/(hi 

. hi+1i  + 

Zytxj^) 

-2hiy(x.  ♦ 

hi*l 

)/(h. 

♦ hl+i* 

“(h?*i  - h 

h2 

ihi+l 

- h?h 

i+l)S'"(xl 

- hi)/6(hi-i-hi+1) 

♦<h*  + 3h. 

hi+l 

* hi*i 

)y  (x.^/6 

-(h?  - h?h 

i+1  ' 

hihi+l)y(xi  + 

hi+l^6*hi*hi+l^ 

now  expand 

y(xt 

‘ hi> 

. y(x.  ♦ 1 

Vl5 . 

M •• 

y (Xj  - hj),  and  y (xi  ♦ hi41)  *n  terms  of  Taylor* s series 
and  then  combine  like  derivatives  of  y,  we  get 
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h.  h 


^xi^  “ ~T£5  ^2hi+l  + 3hihi*l  ' 3hihi+l  " 2h^]y^^(x^) 


h.  h, 


♦ -hfr1  [3h?>i * 2hihi-i  - 7hihi+i " 2hihi+i 

* 3h?ly(6,(x1) 

♦ £’hi*l  * 2hihtl  - ’hihi.l  * ’h?hi.i 

' 2hihi»l  . 

♦ 0(h8)  (It.  3.8) 


We  should  note  that  when  h^  = the  first  and 

third  terms  of  (4.3*8)  are  zero  while  the  second  term  is 
h^y (6) (x^)/24o  and  thus  in  agreement  with  the  Numerov 
method.  However,  when  h^  / hi-*i  the  truncation  error  is 
of  order  5 and,  in  addition,  contains  a ?th  order  term. 

In  Chapter  5 we  describe  a way  to  approximate  T(x^) 
once  we  have  computed  Y,  the  approximate  eigenfunction. 

As  we  did  with  the  other  discretizations,  we  can  now 
introduce  the  matrix  equation 


(A  - XB)  y = T(y) 


(4.3*9) 


where  A is  an  exact  eigenvalue,  y is  an  exact,  discret- 
ized eigenfunction  and 

f(y)  * ( XU^,  t(x2) r(xN)  )*. 

4.3.3  Properties  of  the  Eigenvalues 

We  recall  now,  from  Chapter  2,  that  the  differential 
equation  we  are  approximating  has  an  infinite  number  of 
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real,  distinct  eigenvalues.  Because  of  the  special  sym- 
metry property  of  the  matrix  in  the  Stormer  method,  we 
knew  that  the  matrix  problem  also  had  real,  distinct  eig- 
envalues . 

If  we  expect  our  non-uniform  discretization  to  be 
a good  approximation  then  we  should  expect  its  eigenvalues 
to  be  real  and  distinct.  Certainly  if  we  hope  to  use  the 
Rayleigh  quotient  iteration  to  find  the  eigenvalues  and 
eigenvectors  this  property  is  essential. 

We  are  prepared  now  to  state  the  following* 

THEOREM  4,16*  For  h sufficiently  small  the  N x N 
matrix  problem  (4.3.2)  has  N real,  distinct  eigenvalues. 

Before  proving  this  theorem  it  is  helpful  to  define 
the  following  sequence  of  polynomials. 


P0(*)  - 1 

- (.3 1]L  ( A.  - q.^  - 2)/det(B1) 


Pi  ( ^ ) - ) - 2jdet(Bi_1)pi_1(A) 

- <>i)  * ,12.i-lli;80i<l  - 1i.l>  - ’oi1 

det(Bi_2)pi_2(A ) ‘j/detfB^ 

i- 2, 3 N (4.3.11) 


where  BQ  = 1 and  the  B^,  i-1,2 N,  are  the  leading  prin 

cipal  submatrices  of  B of  order  i. 
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LEMMA  4.17:  The  coefficient  of  X 1 in  p^(X)  is  -*-1 
for  i-0, 1 , . . . ,N . 

Proof  of  Lemma  4.17>  If  we  note  that  detfB^)  - 

then  the  lemma  is  true  for  i=0,l. 

Assume  the  lemma  is  true  for  i=0, 1 j-1 . From 

(4. 3.11), the  coefficient  of  XJ  in  p.(X)  is 

J 

- 1,0jli2.j-ldet(Bj-2)>/det(Bj) 

but  from  relationship  (4.3*6)  for  det(B.)  we  see  that  this 

J 

coefficient  is  equal  to  +1. 

Now  we  note  that  the  N roots  of  pN(X)  •=  0 are  ex- 
actly the  N eigenvalues  of  (4.3.2).  This  is  obvious  if  we 
note  that  the  numerator  of  pN(X), 

^1N^  " qN  ^ " 2 ]det(BN_i)pN_i(X  ) - lf'2,N-l(X  " qN ) ' a2,N-l 

^-30N(X  ’ qN-l^  ~ a0N^det(BN-2)pN-2(^  ' 

is  actually  a 3-term  recurrence  relationship  like  (4.3.6) 
but  for  det(XB  - A)N,  or  more  simply  det(XB  - A). 

Thus  the  roots  of  p^(X)  = 0 are  exactly  the  roots 
of  det(XB  - A)  =0,  which  are  in  turn  the  eigenvalues  of 
(4.3.2).  Thus  the  proof  of  Theorem  4.16  reduces  to  show- 
ing that  Pjj(M  = 0 has  N real,  distinct  roots. 

Proof  of  Theorem  4.16 t The  proof  will  be  by  induc- 
tion. We  will  show  that  p^+^(X)  = 0 has  i-*l  real  roots 
and  that  they  are  separated  by  the  roots  of  p^(X)  - 0. 


i 
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That  is,  between  any  2 consecutive  roots  of  p^^(M  there 
lies  exactly  1 root  of  p^(X). 

Let  i=0.  Obviously  pQ(X)  - 1 a 0 has  no  real  roots. 
Also,  p^U)  = [^^(X  - qx)  - 2}/det(B1)  - 0 has  1 real  root 
at  ( 2 + d n^i ii*  The  separation  requirement  is  satis- 
fied trivially. 

Assume  now  that  p^(X)  = 0 has  i real  roots  s^,  s2, 

...,  s^  separated  by  the  i-1  real  roots  r2 r^  ^ 

of  ~ 0.  i.e.  s-^  < r^  < s2  < r2  < ...  < < s^. 

Then  by  Lemma  4.17  we  know  that 

Pi(A)  = (X  - s1)(X  - s2)  ...  (X  - s^ ) 

Pi_l  ) = (X  - r^^ ) ( X - r2)  ...  (X  - r^) 

Thus 

p^+^(X)  =lLd^^4.^(X  - q^  + ^)  " 2]det(B^)* 

(X  - ...  (X  - Si)  - L^2i^X  ’ qi->-l^  ~ 

a2i-lL'l0,i*l(X  ’ qi)  ‘ ,I0.i*l!de1;(Bi-l)' 

(X  - r1)...((X  - r^^/dettB^) 

(4.3.12) 

If  we  now  let  s.  and  s.4,  be  any  2 consecutive  roots 
pf  p^(X)  - 0,we  see  that  P^+^Sj)  311(1  pi^i^s jt-i^  w111  pe 
of  opposite  sign,  provided  the  factors  »2i^si  " qi+l^  “ a2i 
and  Bq  ^4l(s^  " qi^  ” ao  i+1  are  of  same  sign*  This 
is  adequate  because  (s.  - r^Ms^  - r2)...(sj  - r^  _^)  311(1 
(Sj+l  _ rl^sj+l  ~ 1*2^  • • • ( s j+i  - r^)  are  by  definition 
of  different  signs. 


r 
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The  2 terms  in  question  can  be  assured  to  be  of 
opposite  sign  by  choosing  h sufficiently  small  as  we  did 
in  the  proof  of  Theorem  4.11. 

We  now  know  that  p^^(A)  changes  sign  between  s^ 
and  and  therefore  has  a real  root  in  each  such  inter- 

val. 

We  also  know  from  (4.3.12)  that  since  is  larger 
than  all  the  r^,  wiH  be  negative  when  h has  been 

chosen  sufficiently  small.  However,  Lemma  4.17  implies 
that  p.^(l)  - 00  as  •*  » and  so  Pj^U)  must  have  a real 
root  greater  than  s^ . 

Similarly,  we  know  from  (4.3.12)  that 
has  sign  (-D(-l)1"*  = (-1)1.  However,  again  from  Lemma 
4.17,  P^4^(M  - (-l)*  + ^,0°  as  X - -oo  and  so  P^+j^M  has  a 
real  root  less  than  s^. 

Combining  these  facts  about  roots  we  find  that 
Pi+l(l)  - 0 has  i+l  real  roots  separated  by  the  i roots 
of  p^(* ) = 0. 

This  concludes  the  proof  of  Theorem  4.16. 

Since  the  Numerov  method  can  be  viewed  as  just  a 
special  case  of  our  non-uniform  method,  we  have  immediately 

COROLLARY  4,18 i For  h sufficiently  small, the  matrix 
problem  (4.2.2)  has  N real,  distinct  eigenvalues. 

In  order  to  insure  the  practicality  of  our  non- 


uniform method  we  must  be  able  to  compute  the  upper  bound 
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on  h required  by  Theorem  4.16.  We  must  also  be  able  to 
show  that  such  a restriction  is  not  excessive  for  most 
problems.  Here  we  develop  an  expression  for  such  an  upper 
bound . 

From  Theorem  4.16  we  know  that  h must  be  small 
enough  to  ensure  that 


0o«(s 


2i*sj  " qi*r  ' “2i 


) - a,:  > 0 


(4.3.13) 


for  all  a,  0,  s.,  and  q-t,  and  similarly  we  must  have 

d0,i*l(sj  " qi)  " a0, i+ 1 > °*  (4.3.14) 

We  note  that  conditions  (4.3.13)  and  (4.3.14)  are  suf- 
ficient to  ensure  that  the  sequence  (4.3.11)  is  a Sturm 

sequence.  If  we  let  Q = max  lq(x)|  and  if  we  note  that 

a<x<b 

Sj  > and  thus  let  L = | wQ|  then  by  reasoning  exactly 
like  that  of  Theorem  4.11  we  find  that 


h < 148  / 5 (a  * L)}* 


(4.3.15) 


is  sufficiently  small. 

It  is  important  to  note  that  while  the  condition 
(4.3.16)  is  sufficient  it  is  certainly  not  necessary  to 
ensure  real  eigenvalues.  What  is  necessary  is  that  the 
expressions  (4.3.13)  and  (4.3.14)  have  the  same  sign.  This 
can  occur  in  several  other  ways. 

If  ^2i^sj  * qi*l^  < 0 81,1(1  d o i*l*sj  ~ qi^  < 0 then 
hj  can  be  of  any  size,  and  thus  may  be  very  large.  If 

q^  is  very  small  in  some  subregion  of  La.b]  then  h^  can 


r — 1 
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be  larger  than  (4.3.15)  indicates  since  it  uses  a global 
maximum,  Q. 

We  mention  these  possibilities  here  because  we  find 
in  Chapter  7 that  excellent  results  are  sometimes  obtained 
even  though  (4.3.15)  may  be  violated.  It  is  reasonable 
to  state  then  that  (4.3.15)  is  a sufficient  bound  but  that 
in  practical  situations  h may  be  much  larger. 

4.3.4  Convergence 

As  we  did  with  the  uniform  discretizations,  we  must 
now  consider  the  question!  How  accurately  can  we  approx- 
imate an  eigenvalue  and  eigenfunction  of  the  differential 
equation  with  our  matrix  solution? 

Our  next  result  establishes  convergence  of  the 
matrix  eigenvalue  and  eigenvector  to  the  eigenvalue  and 
eigenfunction  of  the  differential  equation.  To  get  such 
a result  we  appeal  to  the  theory  of  multistep  methods 
for  solving  ordinary  differential  equations. 

We  first  introduce  some  new  notation. 

Let  A be  the  matrix  eigenvalue, 

Y j(  A)  the  jth  component  of  the  matrix  eigenvector, 

Xj  the  differential  equation  eigenvalue, 

Yj(Xj)  the  jth  component  of  a vector  obtained  by 
using  (4.3.1)  as  a multistep  initial  value 
method  as  follows i 

- Let  Yq(X^)  = o and  Y^(X^)  - 1 (arbitrary) 

- Integrate  outward  from  a to  b using  the 


1 


relationship  (4.3.1)  as  a multistep  method 
with  X^  an  exact  eigenvalue  to  get  Yj(X^), 
j=2, 3 N. 

h = max(h. ) , 
i 1 

y(x^)  the  exact  solution  of  the  differential  equa- 
tion at  x^  normalized  such  that  y(x^)  = 1. 

We  now  assume  that  x^,  j > 1,  is  the  first  point 

at  which  the  step  size  changes,  i.e.  h^  - h 2 = ...=  h. 

/ h .+,  . Such  a situation  enables  us  to  use  a well  known 
J x 

theorem  of  multistep  methods  (see  for  example  Gear!  20j, 
p.190.  Theorem  10.8)  to  state  that  both 


Y _i (X i ) = yU^)  ♦ 0(h4) 

and 

Y j (X i)  = y(x  .)  * 0(h4) . 

We  can  now  use  (4.3.1)  to  compute  Yj^(X^)  and  we 

find 


Vi<V  * - 2Yj(xi>  -9ojYj-l(xi> 


-s«Vxi> 

1 / U2j  * 

0 2 j (qi>l 

- X^) 

<< 

X 

u. 

1 

h-» 

oih4)!  - 

2l  y ( x j ) ♦ 

0(h4)] 

-90j(qj-l 

- xi>Lyt*j 

_1)  - 0(hJ 

j 

*>] 

-9lj(qJ  ‘ 

xi)L»txj) 

♦ o(h4)1] 

/ 

U2j  * “20 

(qj*l  * X i 

)} 

- 


but  from  (4.3.7)  we  have 
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y(*Jn)  = l T(xj;i>  - aojyUj.j)  - 2y(Xj) 

' 8o/(xJ-l>  - / 

la2j  " 62jll,j-l  ' V' 

and  so  subtracting  these  equations  we  get 


Yj+l(V  =T(xj)  / la2j  * 02j(Vl  ’ Xi^ 

♦ y(xj+i}  ♦ 0(h4) 

= y(Xj+1)  ♦ 0(h4) 

since  T(x^)  = O(h^)  and  a2j  - 0(1). 

Using  the  numbers  Y.(A.)  and  Y.^.fA.)  as  starting 
values  we  can  again  use  our  recurrence  relationship  (4.3.1) 
as  a multistep  method  to  integrate  to  the  next  step  size 
change  (say  at  xk) . Since  the  errors  in  the  starting 
values  were  0(h  ),  the  theorem  mentioned  above  guarantees 
us  that 

Y^U^)  = y (Xfc)  ¥ 0(h4) • 

This  process  can  be  continued  across  the  entire 
finite  interval  La.bl  so  that  finally 


Vi(V  = y(xN^i}  * 0(h4) 


This  argument  is  only  valid,  however,  if  the  number 
of  step  size  changes  is  finite.  As  h - 0 it  is  theoreti- 
cally possible  that  step  size  changes  could  occur  at  an 
unbounded  number  of  points.  In  practice  we  find  that  large 
intervals  of  uniform  step  sizes  tend  to  develop  in  the 
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adaptive  process  and  thus,  bounding  the  number  of  step  size 


changes  in  theory  will  not  be  restrictive.  We  therefore 
introduce  the  following  definition  and  assume  from  this 
point  on  that  our  method  can  be  so  described. 

DEFINITION » A finite-difference  method  is  uniform 
almost  everywhere  if  there  exists  a constant  K > 0 such 
that  as  h -•  0 the  number  of  points  where  h^  h^ 
remains  less  than  or  equal  to  K. 


We  now  consider  YN<f^(  A)  and  note  thati 

1.  WA)  " 0 

2.  Yn+1(A)  - + (A  - Xi)(dYN+1Ui)/dX) 

♦ o(  A - x^ 2 

* YN*1(V  ♦ * 0(A  " ki)2 

where  zN4.^  * *YN+l^d*  evaluated  at  X^, 

We  must  now  look  at  this  new  sequence  of  numbers 
[Zj}..  Taking  the  partial  derivative  of  (4.3.1)  with  re- 
spect to  X and  evaluating  at  X^  yields 


aOizj-l  + 2zj  4 a2j*j-*-l  * ^0O^(qj-l  " Xi)zj-1  “ yj-l^ 

* *lj[(qj  ‘ Xi)zj  " yj]  + *2j£(qj+l  ’ Xi)zj*l  ’ yj*l]l 
- 0 (4.3.16) 


Writing  (4.3.16)  at  i«l,2 


N we  get  the  matrix  equation 


(A  - XiB)  z - By  (4.3.17) 

In  addition,  since  YQ(A)  ■ 0 and  Yj(A)  ■ 1*  we  know  that 
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*0  " *1  " °* 

Thus  z is  the  numerical  solution  by  our  method  of 
the  differential  equation 

N 

z ■ (q(x)  - X^)  z - y 
z(0)  - z(x1)  » 0 

which  has  ztx.X^)  ■ dyU-^/dX  as  its  solution. 

We  know  from  1.  and  2.  above  that 

0 " *N+1^A^  * YN+l^i^  * ^ * ^i^zN+l  + °^A  “ *i^  ] 

• y(xN+1)  ♦ (Kh^)  + (A  - 

but  y(xN+1)  ■ 0 from  the  boundary  conditions 

.*.  (A  - Vzn-h  * 0(A  " V2  " 0(h4)- 

Because  the  discretization  is  assumed  to  be  uniform 
almost  everywhere , we  know  that  as  h -•  0 the  solution  z, 
of  (4.3.17),  will  approach  the  actual  solution  z,  and  thus 

*N+1  wil1  aPProach  ^N+l’^i^*  Define  f(X)  ■ y(xN4l»M 
as  the  value  of  y(xN+1)  "hen  (4.1.1)  is  integrated  outward 
from  y(0)  ■ 0,  y(x^  ■ 1,  using  4 as  the  eigenvalue.  We 
know  that  when  X - X^,  an  exact  eigenvalue,  ytXj^.X^)  ■ 0. 

If  were  also  zero  then  X^  would  be  a double  root 

of  f(X),  and  X^  would  be  an  eigenvalue  of  multiplicity 
2 which  violates  a property  of  the  S-L  problem.  Thus  if 
we  assume  that  as  h 0,  0(  A - X^)  becomes  negligible 
compared  to  ( A - X^)  we  find  that 
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(A  - a.)  = o(h4). 

Now  since 

VA)  = TjHj)  <•  (A- 

- y(x.)  ♦ 0(h4)  + 0(h4) 

J 

we  also  have 

Y j ( A)  “ y(xj)  + 0(h4). 

We  have  thus  proved  the  followings 

THEOREM  4.19t  If  the  finite-difference  method 
(4.3*1)  is  assumed  to  be  uniform  almost  everywhere,  then 
for  every  eigenvalue,  A,  of  the  differential  equation 
(4.1.1)  there  exists  an  eigenvalue,  A,  of  the  matrix 
problem  (4.3.2)  and  a corresponding  eigenvector,  Y,  such 
that 

( A - 4)  0(h4) 

and  Y.(A)  - y(xi)  0(h4),  i-0,1, . . . ,N+1 . 

We  also  have  as  an  immediate  corollary, 

COROLLARY  4.20s  Theorem  4.19  also  holds  for  the 
matrix  problem  (4.2.2)  on  the  uniform  mesh. 

4.3.5  Error  Estimate 

In  addition  to  the  general  convergence  results 
of  Theorem  4.19,  it  is  possible  to  derive  an  exact  expres 
sion  for  a>4  =4  - A. 


First,  we  need  to  establish  the  fact  that  for  any 
given  A and  for  h sufficiently  small  the  matrix  A - A B 
is  quasi -symmetric . For  this  to  be  true  it  must  be  shown 
that 

Lao,m  * “o,in(qi  " A)]U21  - >»2i(q1*1  - A)] 

can  be  made  > 0 for  all  A and  all  i.  This  is  exactly 
the  condition  required  in  the  proof  of  Theorem  4.16.  Thus 
we  know  that 

h < (48  / «5(U  ♦ | Al))* 
will  suffice. 

This  means  that  there  exists  a diagonal  matrix  P 
such  that  D-1(A  - AB)D  = T is  symmetric  with 


dll  " 1 

dii  = & 26 3 ' ’ ’ 6 i ) ^ i = 2,3, 

* °o j ♦ *oj(qj-i  - A) 

6 j = a2j  4 02j(qj>l  " A)  • 

The  matrix  problem  we  must  solve  is 


(A  - AB)Y  = 0 

or  D”1  (A  - AB)DD_1Y  * 0. 


Letting  D_1Y  - Z,  the  above  becomes 


T7,  0. 


(4.3.1R) 


For  any  pair  (*,y)  of  the  differential  equation 
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(A  - XB)y  * r(y) 

or  D”1 ( A - XB)DD’^y  - D*1  f (y) . 

Multiplying  by  Zt,  we  get 

zVX(A  - X B)DD_1y  = zV1  f(y).  (4.3.19) 

Letting  y - Y ♦ Ay  and  X - A + AX , and  substitut- 
ing into  (4.3.19)  we  get 

ZtD’1tA  - AB  - axb } D ( Z + D_1Ay)  = Z^-1  T(y)  . 

Multiplying  out,  the  expression  becomes 

ZtD’1(A  - AB)DZ  ■»  Z^'-^A  - AB)DD_1Ay  - ZtD"1AXBDD'1Y 
- Z^D_1aX BDD_1Ay  = Z^'1  f (y)  . (4.3.20) 

The  first  term  of  (4.3.20)  is  zero  directly  from  (4.3.18). 
The  second  term  is  zero  because 

ZtD'1(A  - ABjDD^Ay  = ZtT(D“1Z.y), 

but  TZ  = 0 so 

(TZ)t  = zV  =■  ZtT  = 0. 

Thus  from  (4.3.20)  comes  the  identity 

AX  = - zV1  f(y)  / ZtD"1By.  (4.3.21) 

We  have  assumed  here  that  Z^’^By  / 0.  This 
assumption  is  assured  if  the  following  two  conditions  are 


r 
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met. 

First,  the  term  YtD”1D"^By  must  not  be  inherently 
zero  due  to  orthogonality.  Since  the  eigenvector  can  be 
normalized  in  any  fashion,  we  choose  to  require  that 
YtD~*D“1BY  = 1.  From  Theorem  4.19  we  now  see  that 

YtD"1D"1By  = YtD'1D"1B{Y  ♦ 0(h4)} 

= 1 + 0(h4) 

and  therefore  not  equal  to  zero. 

Second,  the  mesh  points  lx.}  cannot  all  be  at  nodes 

J 

of  the  eigenfunction  y(x)  or  the  term  in  question  would 
again  be  zero.  This  only  requires  that  h be  small  enough 
to  ensure  that  when  finding  the  kth  eigenvalue  there  are 
more  than  k-1  mesh  points.  This  requirement  is  necessary 
in  any  event  because  to  find  the  kth  eigenvalue  with  any 
accuracy  our  matrix  problem  must  at  least  be  of  size  k. 

If  we  make  the  additional  assumption  that  the  fourth 
term  of  (4.3.20),  which  contains  the  product  AXAy,  is 
negligible  compared  to  the  other  terms,  we  get 

AX  * - f(y)  / ZtD“1BY 

= - YtD’1D"1  ?(y)  / YtD‘1D'1BY.  (4.3.22) 

Since  Y is  known  while  y is  not,  (4.3.22)  provides  an 
approximation  for  the  error  in  the  matrix  eigenvalue  which 
can  be  computed.  This  estimate  is  used  in  Chapter  5 as 
part  of  our  adaptive  method. 
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4.3*6  Solution  of  the  Matrix  Problem 

Again,  as  with  the  two  uniform  methods,  we  use 
Rayleigh  quotient  iteration  to  solve  the  matrix  problem 
(4.3.2).  Since  A and  B are  no  longer  symmetric  we  can 
expect  only  linear  convergence,  according  to  Ostrowski, 
with  Algorithm  4.5*  In  order  to  get  cubic  convergence 
with  general  A and  B we  must  use  the  generalized  Rayleigh 
quotient  in  which 

A - v*Au  / v*Bu 

The  u and  v are  right  and  left  eigenvectors  respective- 
ly. Ostrowski  goes  on,  however,  to  show  that  if  the  prob- 
lem has  real  eigenvalues  with  linear  elementary  divisors 
then  the  simple  Rayleigh  quotient  iteration  (Alg  4.5) 
produces  quadratic  convergence.  He  argues  that  since 
finding  both  left  and  right  eigenvectors  takes  approx- 
imately twice  as  much  work,  we  can  perform  2 iterations  of 
the  simple  scheme  for  every  generalized  iteration.  Thus 
the  simple  method  gives  effectively  fourth  order  con- 
vergence when  compared  with  the  generalized  method. 

Since  by  Theorem  4.16  we  know  that  our  problem  has 
real,  distinct  eigenvalues  when  h becomes  small  enough, 
and  therefore  has  linear  elementary  divisors,  we  again  use 
Algorithm  4.5  to  solve  the  matrix  problem  (4.3.2).  Al- 
though quadratic  convergence  is  expected,  actual  numer- 
ical tests  show  that  cubic  convergence  is  again  obtained 
during  at  least  some  stage  of  the  iteration  process. 


70 

In  tables  4.3  and  4.4  we  present  rate  of  convergence 


results  for  the  following  problem. 

y ♦ (X  - x|x|)y  = 0 

y(  -1)  = y( 1)  - 0 

Table  4.3  is  for  XQ  » 2.462*58070  and  Table  4.4  for  A ^ » 
22.2077749. 

In  each  table  A^  was  computed  by  varying  sizes  of 
non-uniform  meshes.  For  each  matrix  size  N,  the  Rayleigh 
quotient  iteration  was  performed  until 


^i.k  “ A-i 


i,k*-l 


I < 10 


-13 


H was  then  computed  from  the  equation 


i .k 


i . k - 1 


I \ht 


i.k-1 


- Ai.x' 


We  see  from  Tables  4.3  and  4.4  that  cubic  conver  - 
gence  or  better  is  obtained  for  each  matrix  at  some  point 
in  the  iteration.  This  can  be  explained  by  the  close 
association  between  our  matrix  problem  and  the  differen- 
tial equation. 


■ 


J 


■■ 


Table  4.3 

Convergence  Results  for 


1 = 28 

= 30 

N = 32 

N 

- 34 

: & 

k 

3 

k 3 

k 

3 

2 

2.11 

2 

2.27 

2 

2.40 

2 

2.46 

3 

3.67 

3 

3-55 

3 

3-37 

3 

3.26 

4 

2.39 

4 

2.33 

4 

2.29 

4 

2.28 

5 

1.65 

5 

1.70 

5 

1.77 

5 

1.74 

Table  4.4 

Convergence 

Results 

for  * 

2 

N 

= 8 

N 

= 16 

N 

- 32 

k 

3 

k 

3 

k 

3 

2 

1.79 

2 

8.23 

2 

11.1 

3 

, 3.18 

3 

2.48 

3 

1.91 

4 

1.49 

N 

* 49 

N 

* 61 

N 

= 73 

k 

3 

k 

3 

k 

3 

— 

2 2.28 

3 3.?0 

4 2.29 


2 2.27 

3 3.51 

4 2.26 
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Our  problem  can  be  written  in  operator  notation  as 
Ly-^y  (4.3.2?) 

where  L ir  the  Sturm-Liouville  operator 
!jU  - -(p(x)u* ) «•  q(x)  u J 

p,  p , q,  and  r real -valued  and  continuous,  and 
J'(x)  > 0,  r(x)  > 0 for  a < x < b. 

Thus  LiU  - H where 

b ? 

H - i.u(x)  i J |u(x)rrlx)  dx  < »}  and 

a 

2 

l)  - (u(x)  s u f C (a  < x < b)  which  satisfy  the 
boundary  conditions) 


and  where 

b 

lu,v)  = J u(x)  v(x)  r(x)  dx. 
a 

We  can  now  define  Rayleigh  quotient  iteration  for  (4.?. 23). 


1L  - A ^r(x)  - r(x)y^  (4.3.24) 

x(k*i)  . j-f<*i)Ll(it’i)<|X  / (2(kn)tl(^i)) 


.(k+1) 


,(k  ► !) 


where  v is  some  normalization  factor.  All  integrals  are 
assumed  to  be  over  (a,b). 

Since  yv  and  zv  are  in  U,  and  L is  the  S-I 

operator,  it  is  well  known  that  each  has  the  expansion 


y[K) (x)  = E a.  y. (x) 


z(k+1)  (x)  = £ b.  y.  (x) 

i=l  1 1 
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(4.3.25) 

(4.3.26) 


where  the  ty^(x))  are  the  eigenfunctions  of  L and  where 
a^  = (y^k\y^)  and  b^  = (z^k+3^,y^).  All  £ are  assumed 
to  be  over  1 to  *. 

Since  the  normalization  factor  v can  be  chosen 

( k) 

arbitrarily,  we  assume  that  y has  been  normalized  such 


that 


r(k)  s v . ♦ £ ' £ , y. 


i *i 


where  £ denotes  and  = €*d^  - 0(0  . (i.e.  each 

i/j 

coefficient  is  assumed  to  be  a multiple  of  some  small  num- 


ber ( . ) 


We  wish  to  show  that  after  one  iteration  of  this 


process 


y(k+1)  - y.  ♦ e'  o ( € 3 ) Yi. 

Substituting  (4.3.25)  and  (4.3.26)  into  (4.3.24)  we  get 


Thus, 


1L  - Xmr}  Eb^  “ r Ea^ 


LL^biyi  ] - K K;r  Eb^  •-  r Ea^ 


Now  because  the  operator  L is  a bounded,  linear  operator 
we  can  write 
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EbiLyi  - *(k)  Sbiryi  = Ea^y^ 

But  Ly^  - k^ry^,  from  which  it  follows  that 
Lbi  - A-(k))ryi  = Ea^y.. 

and  since  the  eigenfunction  expansion  is  unique  we  have 
b^  = a ^ / (A^  ~ ^)*  i-1 , 2, . . . 


Thus 


z - Eb.y 


Fi 

j ' Uj 


- y,  / U,  - *(k))  • t’yil  ^ / (^i  - *(k)>J 


Then  when  z^k*^  is  normalized  to  yield  y^^we  find 


y(k*1}  - yj  - 


(xj  - l{k))*'yiLti  / - *(k))  1 

(4.3.27) 


Now 


A(k)  ..  jy(k)Ly(k)dx  / (y(k,.y(k)). 


( k) 

Substituting  the  series  for  y we  find  that 

A(k)  £ilaiaj4jj'y1ryjdx  / lia^J  y^y^dx 

Because  of  the  orthogonality  of  the  eigenfunctions,  i.e. 

(yi #y  .i ) = 6 i i • we  havf? 

A(k)  i:a?  A . / Ea? 

• *=’•?  *-4l  / LI  • *•«?! 

- 2..  * . 2 v l / ii  . r 2r  1 j2 


* U . d£  Xj  J / 11  ♦ df]. 


(4.3.2ft) 
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Now  because  the  series  (4.3*25 ) is  convergent,  the  series 
in  both  numerator  and  denominator  of  (4.3.28)  converge 
and  thus 


x(k)  = [x  . ♦ o(e^)  ] / LI  * o ( t ^ ) ] 

J 

r X . 4 0(c2)  . 

J 


-2, 


(4.3*29) 


Thus  substituting  into  (4.3.27)  we  find  that 
y(krl)  - y 


-(ka)  - » ♦ F.’y.  0(t3) , 


.1 

showing  that  cubic  convergence  is  obtained. 

We  note  in  Tables  4.3  and  4.4  that  cubic  convergence 
is  indeed  obtained.  We  also  note,  however,  that  in  every 
column  the  convergence  rate  attains  the  cubic  level  and 
then  immediately  drops  off  toward  unity. 

( k ) 

Suppose  that  instead  of  exact  solutions  y we 
can  only  find  approximations 

~(k)  . h^e(x)  --  y^  * O(h^) 


where  e U,  i.e. 


Then 


i ( x ) - Ec^U) . 


.00  . * h>m  vu|  * A** 


(4.3.30) 


J [>  ( a^  h ci)yi  (ai  ♦ h c^ry.^  dx 
/L ^ ( a^  «•  h^c^)y^ JL^(a^  + h c^)ry^|dx 


1 


' 

[ LI 
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» . 

SEta^  ♦ h c^)(a^  ♦ h c.i)^.ii,yiry.1  dx 

EEta^  + h c^Haj  ♦ h c^/y^ry^  dx 

E(a,  ♦ h^c*)2*,  dx 

3 ■ * ■ — - ■ * 

E(a^  + h4ci)2  dx 

Ea2^  ♦ 2h4EaiciXi  ♦ h8Ec2Xi 
Ea2  ♦ 2h  ta^c^  ♦ h Ec2 

I ^ 

and  thus  since  (4.3*25)  and  (4.3.30)  converge  we  have 

k(k)  . Vi  * olf2>  * 

1 ♦ 0(€2)  4 0(h4) 

» Xj  ♦ 0(c2)  •*•  0(hu). 

Substituting  this  into  (4.3.27)  we  see  that 

y(R+1)  " yj  4 2'yi[o(€2)  ♦ 0(h4)]fi  / (^  - X(k)) 
Thus  for  any  fixed  h we  find  that  while 
0(h4)  < < 0(€2) 

the  convergence  is  cubic  as  before,  however,  after  several 
iterations  when 

0(h4)  > 0(€2) 

A j. 

the  rate  becomes  simply  linear  because  LoU*)  ♦ 0(h  )]€^ 
is  0( € ) . 

This  phenomenon  is  borne  out  by  the  results  of 
Tables  4.3  and  4.4. 

I I ri 

E,:  r T 


. 


4.3.7  Deferred  Correction 
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With  the  use  of  the  estimate  (4.3.22)  for  aA  and 
with  an  approximation  T ( Y ) to  t(y)  we  can  again  apply 
the  deferred  correction  to  improve  our  final  matrix 
solution. 

The  algorithm  we  use  is  as  follows. 

ALGORITHM  4.21i  Let  Y and  A be  the  solution  to 
the  matrix  problem  (4.3.2). 

1.  Estimate  ?(y)  by  computing  T(Y) 

2.  Estimate  aA  by  computing 

aT  = - ZtD"1T(Y)  / ZtD"1BY 

3 . Le  t S’  ~ A + aT 

4.  Solve  (A  - ^B)y  = T(Y)  for  y 

If  aA  were  found  exactly  in  Step  2 using  expression 
(4.3.21)  rather  than  (4.3.22)  then  we  could  compute 
A = A ♦ aA  and  y exactly.  However,  we  must  use  instead 
the  expression 

AT  = -YtD"1D’1f(Y)  / YtD"1D’1BY 

The  errors  in  this  expression  occur  in  two  places. 
First,  concerning  the  second  Y in  the  denominator,  we 
know  from  Theorem  4.19  that  Y = y 0(h  ).  Second,  there 
is  an  error  when  T(Y)  is  used  to  approximate  T(y). 

We  have  not  yet  introduced  an  approximation  T(Y) 
to  t(y) , choosing  rather  to  present  it  along  with  the 
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practical  details  of  the  adaptive  method  of  Chapter  5. 

We  assume  for  now  that  such  an  approximation  is  available 
and  that 

T(Y)  = t(y)  ♦ 0(h°) 

recalling  that  t(y)  = O(h^) . 

With  these  two  sources  of  error 

£X  = -YtD~1D”1[  f(y)  + 0(h8)]/YtD"1D"1BLy  + 0(h4)] 

At  this  point  we  note  that  as  h - 0 the  number  of 
mesh  points  N - as  l/h.  Thus  YtD"^D"1BY  = 0(h)  and 
YtD“1D”1[o(h8) ] = 0(h7)  and  so 

sx  , , 4 Sibil 

YtD'1D”1By  ♦ 0(h5)  0(h) 

= AX  [\  + 0(h4)}  ♦ 0(h6) 

= Al  + O(h^)  since  ^ = 0(h4). 

/.  A ♦ aT  = A + aX  ■*  0(h6) 

= A + 0(h6) 

. k 

Thus  we  have  shown  that  while  A - A ♦ 0(h  ) one 
correction  step  gives  us  sixth  order  accuracy  as  before. 

The  non-uniform  finite-difference  method  developed 
in  this  section  becomes  the  basis  of  our  adaptive  Sturm- 
Liouville  solver  of  Section  5*2.  Numerical  results  from 
the  method  are  presented  after  the  adaptive  algorithm  is 
discussed . 
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The  contributions  of  this  chapter  have  been  to 
provide  us  with  i)  a finite-difference  method  which 
uses  the  non-uniform  mesh  we  expect  from  an  adaptive 
method,  ii)  a cubically  convergent  technique  for  solving 
the  resulting  matrix  problem,  and  iii)  a process  for 
correcting  the  matrix  solution  which  increases  the  accur- 


acy from  O(h^)  to  O(h^) . 


m 


Chapter  5 
AUTOMATIC  METHODS 

5.1  Introduction 

Davis  and  Rabinowitzi.12 ] , in  their  book,  define 
four  categories  of  numerical  integration  methods.  These 
definitions  also  apply  to  finite-difference  methods  for 
solving  the  Sturm-Liouville  problem.  They  are: 

Adaptive  - A finite-difference  solution  method  is 
adaptive  ::.f  the  points  of  the  mesh  are  chosen  to  depend  on 
the  nature  of  the  solution,  (y.M. 

Non -adaptive  - A method  is  non-adaptive  if  the  mesh 
is  chosen  according  to  a fixed  scheme  independent  of  the 
solution . 

I terative  - A method  is  iterative  if  successive 
approximations  to  the  solution  are  made  until  agreement 
within  a prescribed  tolerance  is  achieved. 

Non-iterative  - A method  is  non-iter-.tive  if  infor- 
mation taken  from  the  first  approximation  is  used  to  gener- 
ate a second,  which  is  then  taken  as  the  fina1  result. 

With  these  definitions  we  are  able  to  categorize 
the  methods  of  Chapter  4 as  non-adaptive , non-iterative 
methods.  The  mesh  is  chosen  before  using  the  algorithms 
and  then  the  deferred  correction  is  performed  once  to  get 
a final  result. 
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In  this  chapter  we  present  two  iterative  methods,  a 
non-adaptive , iterative  method  in  Section  5*2  and  an  adap- 
tive, iterative  method  in  Section  5.3* 

Both  methods  require  accurate  initial  estimates  to 
ensure  convergence,  and  neither  method  can  be  used  to  solve 
singular  S-L  problems.  Such  a method  we  define  as  non- 
automatic . We  say  that  a method  is  automatic  if,  when 
supplied  with  only  the  differential  equation,  boundary 
conditions,  and  a desired  tolerance  level,  it  can  return  an 
approximate  solution,  (Y.A),  accurate  to  within  the  spec- 
ified tolerance  on  a. 

In  Chapters  6 and  7 we  address  these  shortcomings. 
The  final  result  is  an  automatic,  adaptive  finite-differ- 
ence method  for  solving  the  S-L  problem. 

6.2  A Non-adaptive.  Iterative  Method 

In  this  section  we  present  a non-automatic , non- 
adaptive,  iterative  method. 

6.2.1  Background 

The  algorithm  to  be  introduced  is  motivated  by  the 
work  of  Lontini  and  Pereyral26,2?J.  In  these  papers,  a 
variable  order  finite-difference  method  is  developed  for 
solving  nonlinear  multipoint  boundary  value  problems.  Eig- 
envalue problems  are  not  considered. 

The  general  problem  solved  in  that  work  is  the  non- 
linear, first  order  system 

y (t)  - f(t,.y(t)  ) “ 0,  a < t < b 
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subject  to  the  multipoint  boundary  conditions 


g ( y ( t^ Jtyttg)*  • • • > y ( t^ ))  ■ o , a < t^  t^  £ b . 


The  method  is  based  on  the  technique  of  deferred 
corrections.  It  is  assumed  that  the  following  asymptotic 
expansion  for  the  local  truncation  error  is  available. 


W 


K 

E 

k-1 


aky  <2k,2>(*.)h2k  * 0(h2K<2> 


(5.2.1) 


The  first  m terms  of  (5.2.1)  can  then  be  approximated  by 
a linear  combination  of  function  values  at  various  points 
as  such  1 


T 

m 


m 

(x.)  - - I 
1 k=l 
2m*  2* 
= E 
s=l 


aky  (2k42>(x1)h2k*2 
q 

wsy  (xj  * agh)  t 0(hq) 
(x.))  * 0(hq) 


(5.2.2) 


where  the  a are  integers.  The  numerical  approximation  is 

5 

found  by  discretizing  the  differential  equation  and  solving 
the  resulting  system. 

If  Fh(Y)  = 0 is  defined  as  the  nonlinear  discretized 
problem  and  Y « (Y^,Y2 as  the  discre-tized  solu- 

tion, we  find  that  Pereyra  in  L32]  introduces  the  follow- 
ing iterated  deferred  correction  algorithm. 


ALGORITHM  5. It 

1.  Let  Y^  be  an  0(h21t' ^discrete  solution. 

2.  Compute  h-2Sk+^ ( Y ^ ) , an  0(h2**2)  approximation 
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to  the  first  k+1  terms  of  the  local  truncation 
error . 

3.  Solve  Fh(Y)  - h‘2Sk+1(Y(k))  for  Y(k4l) . 

It  is  obvious  that  only  a limited  number  of  iter- 
ations of  Algorithm  5.1  can  be  performed  on  a mesh  of  size 
N,  because  an  appropriate  number  of  mesh  points  must  be 
available  surrounding  each  point  to  form  Sk-i>1(Yv 
Pereyra  also  notes  that  for  a given  problem  and  mesh  size, 
only  a limited  number  of  corrections  will  produce  improve- 
ments to  the  solution.  As  Pereyra  puts  it,  this  is  due  to, 
"the  growth  of  high  order  differences". 

Equipped  then  with  an  estimate  to  the  error 

( k-1 ) — 

e^_^  - Y - y , Lentini  and  Pereyra  introduce  the  fol- 

lowing iterative  algorithm. 

ALGORITHM  S.2« 

1.  Let  k 0. 

2.  Is  N > 4 ? 

NOs  Refine  mesh  and  GO  TO  1. 

YES*  Then  S^(Y^0^)  can  be  computed,  so  solve 
Fh(Y)  - 0 for  Y(0),  compute  S^Y10*), 
and  compute 

3.  If  II  aq  ||  < TOL  (the  user  provided  tolerance) 
then  STOP. 


Correction  Loop: 
4.  Let  k k'l . 
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5.  If  N < 2k+2  then  refine  the  mesh  and  GO  TO  1. 
(not  enough  points  to  compute  SR) 

6.  Solve  F.  (Y)  = S^(Y(k“1))  for  Y(k). 

n k 

7.  Compute  (Y  ^ k^ ) . 

8.  Compute  ar. 

9.  If  II  A ||  < TOL  then  STOP. 

10.  If  II  ar  ||  < C II  ar_1  II,  (0  < C < 1).  then  GO 
TO  4 else  refine  mesh  and  GO  TO  1.  (This  en- 
sures that  the  errors  are  decreasing,  i.e. 
it  checks  for  growth  of  high  order  differences) 

Pereyra  admits  that  there  are  some  theoretical  dif- 
ficulties in  obtaining  the  expansions  necessary  to  justify 
the  method.  Different  differentiation  formulas  must  be 
used  at  different  points.  This  problem  is  at  least  delayed 
by  replacing  higher  derivatives  of  y by  derivatives  of 
f(x,y)  two  orders  lower. 

A further  difficulty  we  encounter  when  trying  to  use 

Algorithm  5*2  to  solve  our  problem  is  that  Pereyra  only 

considers  uniform  spacing  when  computing  the  coefficients 

w of  (5*2.2).  While  this  is  sufficient  for  the  present 
s 

non-adaptive  algorithm,  it  is  not  for  our  ultimate  adaptive 
method.  Lentini  and  Pereyra  in  l.  26  | and  ( .27  1 do  consider 
a non-uniform  mesh,  however,  there  again  only  first  order 
systems  are  discussed. 

We  should  also  recall  that  our  objective  is  not. 


as  was  Lentini  and  Pereyra’ s,  to  develop  a variable  order 
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method.  We  rather  want  to  keep  the  order  of  the  method 
constant  and  increase  accuracy  by  refining  the  mesh  in  the 
proper  way. 

With  these  difficulties  and  considerations  in  mind, 
we  choose  to  introduce  a simpler  non-adaptive  algorithm 
which  can  be  easily  converted  to  an  adaptive  method. 

5.2.2  Truncation  Error 

Our  first  task  is  to  develop  an  approximation  to  the 
local  truncation  error.  The  basis  for  our  non-adaptive 
algorithm  is  the  Numerov  method  of  Section  4.2.  The  non- 
uniform  finite-difference  method  of  Section  4.3  becomes 
the  basis  for  our  adaptive  algorithm  of  the  next  section. 
Since  the  Numerov  method  is  simply  a special  case  of  this 
non-uniform  method  we  develop  here  the  approximation  neces- 
sary for  the  more  general  problem. 

We  recall  from  Section  4.3  (equation  4.3.8)  that 
the  truncation  error  can  be  written  as 

t(xt)  = y(5) (xi)hihi+1ai/l80 

♦ y(6)(x.)hihi+1b./720 
4 y(7)(xi)h.hi+1c./5040 

♦ 0(h8)  (5.2.3) 

•i  = 2hiU  * 3hihi+i  - 3h2h.+1  - 2h3 

bi  = 3hi*l  * 2hihi*l  ' 7hihi+l  > 2hiVl  + 3hi 


where 
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c . = 


5hi+l 


2hihi*l 


9hihin 


9hi»lU 


■2hihiU  - 5hr 


We  know  from  the  derivation  of  the  non-uniform 
method  that  it  is  exact  for  y(x)  = xp,  p=0,l,2,3,4,  and 
thus  T’(Xj)  = 0 when  y(x)  is  of  this  form.  Now  let  us 
compute  t^  = T(x^)  when  y(x)  = xp,  p~5.6,7.  (or  equiv- 
alently y (x)  = p(p-l)xp'2,  p=5,6,7) 


ti5  ""  5!hih 

because  y ^ (x)  = 5s 
derivatives  are  zero 


lai/18°  - 3h.hi+1a. 


when  y(x)  = x^  , while  all  higher 
. Similarly, 


ti6  » 6!xihihi+1ai/180  + 6!  hjh^^/720 


9xihihi+iai  - 
M v ^ 


and  ti?  = 7!xfhihi41ai/360  « 7!  xihihi+lbi/72° 

f?!hihinci/504° 

* 14xihihinai  ' 7xihihi.ibi  ' hi*Vi°i 


We  now  assume  that  the  approximate  solution,  A and 

M 

Yit  i=0,l, . . . ,N+1,  has  been  computed.  Prom  it  Y^  can  be 

N M 

computed  by  Y^  = (q(x^)  - A )Y^.  Letting  f^  = Y^ , we  form 
the  following  interpolating  polynomial  on  the  points 
tXj,  j=i-3.  i-2.  i-1,  i.  i+1,  i+2}i 


gi(x)  = fi  + (x  " xi)fL^i_1,xiJ  ♦ (x  - xi_1)(x  - x.)- 
^-xi-l'xi ,xi+l 1 * (x  - xi_1)(x  - x±) (x  - 
f^xi-2,xi-l,xi,xi+l^  ■*  (x  * xi_2)...(x  - xi+1)« 
ftxi_2'  * * *xi+2  ^ ^x  " xi-2^***^x  " xi+2^* 

ft-xi_3*  • • ,xi+2^ 

= Aj  + B jX  + C^x2  ♦ DjX3  E^4  ♦ F^5 

where  f£x ^ ,x , . . .x is  the  kth  divided  difference  of 
•• 

Y at  the  points  x^,  xj+i»»»»»  x^^. 

Solving  for  D^,  E^,  and  F^,  we  find 

Di  = f^xi-2*  * * *xi+l  ^ “ (xi-2  + xi+l)* 

f[xi_2, . . .xi+2]  + sif[xi_3, . . ,xi  + 2J 

where  s.  = (xi.2xi-l+xi.2xi"xi-2xi>l+xi-2xi*2"xi-lxi+ 
xi-lxi4‘l^xi-lxi*  2fxixi+l*xixi+  2+xi*  lxi+2) 

Ei  = xi-2’ * * *xi+2-^  " ^xi-2  * * * * *xi+2^  * 

fLxi-3  • • • ,xi+2^ 

Fi  = ftxi-3*  • • ,xi+2-^ 

N 

Regarding  g.(x)  as  an  approximation  to  y (x)  on  the 
points  tx^_3 xi<>2J , we  can  now  approximate  the  trun- 

cation error  TU^)  by  substituting  g-^x)  into  (5.2.3) 

H 

for  y (x).  This  yields  the  following  approximation  zo 
t (xi) . 

t<*i>  * Ti(xl>  * iKs  1 * 12*17 


( *>.  ?.4) 


n 
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For  the  several  points  at  the  ends  of  the  interval 
where  . . .x^  g I cannot  be  evaluated,  we  choose  to 

approximate  T(x^)  by  extrapolating  the  last  computable 
T^.  For  example, 

D_  E F 

T(x2)  =»  T^(x2)  ~ 20^25  * 30t26  ’ 1*2*27  (5*2.5) 

Just  how  good  an  approximation  we  have  is  established  in 
the  next  theorem. 

The  divided  differences  which  appear  in  the  coef- 
ficients D^,  E^,  and  F^  of  (5.2.4)  are  divided  differences 

••  m m m h 

of  Y , not  y . Theorem  4.19  shows  that  = y (x^)  ■»  0(h  ) 
If  the  errors  in  the  Y^  were  just  random,  we  would  find 
that  we  lose  an  order  of  h for  each  divided  difference 
taken.  For  example, 

,r.  - , _ Yi  - h-1  - y"(xl)  - * Olh^i 

nxi-l,xiJ  " hi  " hi  hi 

= y Lxi.xi_1l  + o(h^) . 

M 

But,  if  the  errors  are  assumed  to  be  smooth,  i.e.  Y^  - 

M li 

y (x^)  + e(x^)h  where  e(x)  has  at  least  as  many  deriva- 

m h 

tives  as  y , then  all  divided  differences  are  still  0(h  ). 


flxi  rx 


_ Yj  - y1,1  _ y ty  - ,v  (xi_1) 

i ’ h. 

4 eUj)  - e(xi,1) 


hi 


= y + e (x^n*  + o(h^) 

- y l xi_i»xi  1 4 o(h4) . 


Uith  this  assumption  we  get  the  following* 


L 


' 
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THEOREM  5 . 3 « Let  y(x)  be  the  solution  to  (4.1.1) 
and  be  8 times  differentiable  on  [a,b]  and  let  = 

H II 

y (x^)  ♦ e(x^)h  , where  e(x)  is  6 times  differentiable. 

Then  T^(x^)  " T(x^)  = 0(h8),  where  h = max(h^). 

Proof > We  let  f^tx^)  be  the  quantity  (5*2.4)  com- 
puted with  y L...J  rather  than  fL  . * . .] . From  the  above 
discussion  and  the  assumption  on  e(x),  we  know  that  = 

♦ O(h^),  and  similarly  for  E^  and  . (the  coefficients 
of  Ti(xi)  ) 

Thus,  since  t^  = O(h^),  t^^  = 0(h8),  and  tj^  = 

n 

0(h),  we  see  that 

T.(xi)  = T^(x^)  + 0(h9).  (5.2.6) 

We  now  find  the  error  in  T^(x^).  Our  procedure  for 
finding  T(x^)  is  to  fit  an  interpolating  polynomial  to 

N _ N 

the  points  [(x^y  ( xi ) } , namely  g (g^,  using  v (,...! 
rather  than  fL...])  From  Conte  and  deBoorlloJ  we  know  that 

M — H _ 

y (x)  •-  g^(x)  ♦ y Lxi_^j, . . ,xi42,x  Jp^(x) 

= ii(x)  ♦ y^8) ( 0)P6(x)/6! 
i+2 

where  P,(x)  « t»  (x  - x.)  and  0 € Lx.  ,.x..0], 

6 j=i.3  ) i-3  1+2 

— •• 

Thus  when  g^(x)  is  substituted  into  (5*2.3)  for  y (x), 

we  must  compute  its  derivatives,  for  example 


! 


fT 

li 

T- 

1: 


(dxy"'-xi-3’ 


X 


.Xif 2,x  ]P6(x)  ♦ 
xl42.x 1p^(x) . 


i-3'  * * * 


II 


Since  T^(x)  will  only  be  evaluated  at  a mesh  point. 


Pg(x^)  “ 0 and  80 

g[(x)  - y(?)(x)  - y(8)(91)P^(x)/8! 

At  any  mesh  point,  x^  € Lxj^.x^^] 

. i*2  , 

P,(x.)  - it  (x.  - x.  ) ■ 0(h5) 

6 J k-i-3  J k 


and  so 


g!  (x  .)  * y^  (x  .)  ♦ 0(h5)  . 

* J i' 

Continuing  to  take  derivatives  of  g^(x)  we  find  that 

i[k)(xj)  - y(k+2)(xj)  ♦ 0(h6-k).  k-1,2 5 

j-i-3. . . . ,i+2 

i.e.each  derivative  of  g^  loses  an  order  of  h in  accuracy. 


So  we  have 


W " £i{3)(xi>hihi4iai/180  ♦ 
ii(4)(x.)hihi+1bi/720  ♦ 

g{(  5)  (xi)hihi^.1ci/5040 

■ ty^^(Xj_)  ♦ CXh^nhjhj^aj/180  ♦ 

[y(6)(xi)  ♦ 0(h2)]hihi^1bi/720  + 

[y(7)(xi)  ♦ 0(h)'Jhihi+1c1/5040 

- t (xL)  ♦ 0(h8). 
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Combining  with  (5-2.6)  we  have 


T-(x.)  - t(x.)  ♦ 0(h°) 


and  the  theorem  is  proved. 

5.2.3  The  Non-adaptive.  Iterative  Algorithm 

We  now  present  a non-adaptive,  iterative  algorithm. 
We  recall  the  following  from  Section  4.2. 

1.  The  discretization  (4.2.1)  (Numerov  method) 
leading  to  the  matrix  problem  (4.2.2): 


AY  - A BY  0 . 


2.  The  truncation  error: 


T(x^)  - 2%q  h6y^(0). 


3.  The  error  estimate  (4.2.6): 


aA  « - Y*  t(.v)  / YtBY. 


[ 


r 


n 


n 


4.  Algorithm  4.6,  used  to  correct  the  matrix  eig- 
envalue and  eigenvector. 

We  also  note  that  when  h^  - h^+^  for  all  i , equation 
(5-2.5)  reduces  to 

W ' 255  <*1-2  - W"l-1  * 6*I  - "Vi  * *1.2) 

- h!  v(6)m 
250  y * 

This  is  easily  shown  as  follows. 


When  hj^  - hi+1  ■ h,  ti$  - 0,  ti6  - 3h  , and 


ti?  - 21h  so 

Ti(xi*  " TO  h6  * 2^xih 

h6 

* 10  ^xi-2’  • * *xi+2-^  “ ^xi-2  'f,,,+xi+2^* 

flxi_3*  * • ,xi*2-l  * 5xifLxi_3*  • • *xi^ 

but  since  h^  ■ hj^  ■ h for  all  i» 

(xi-2  * *i-i  + xi  * xm  * xi+z’  ' 5xf 
6 

Ti(xi)  * Yo  f^xi-2*  * * ,xi+2^ 

- 15  ^ 

where  A1  is  the  ith  forward  difference 

- 2^  C*"i-2  - <-l  * 6YI  ' 4YI-1  + YW>- 

With  these  preliminary  remarks  we  are  now  able  to 
present  the  following  non-adapti ve,  iterative  algorithm. 

ALGORITHM  5.4i 

1.  Set  up  a basic  uniform  mesh  wQ  on  [a,b].  (We 
use  N«8  in  most  cases.) 

2.  Solve  the  matrix  problem  (4.2.2)  using  Algo- 
rithm 4.5* 

3.  Compute  error  estimate 

a\  . _ Y*  TO)  / Y 1 HY 


and  make  one  correction  using  Algorithm  4.6 


r 


I 


n 
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4.  If  |aT|  < TOL,  STOP. 

5-  define  mesht 
i)  N = 2N 

ii)  Interpolate  values  of  Y^  at  new  mesh 
points  to  get  starting  right-hand  side  for  next 
iteration. 

6.  GO  TO  2. 

Note » Step  5ii)  is  not  necessary  but  it  speeds  up  con- 
vergence of  Algorithm  4.5.  After  the  first  iteration 
(i.e.  N > 8),  Algorithm  4.5  is  changed  slightly.  Rather 
than  letting  BYQ  = (l,l,...,l)t,  the  algorithm  computes 
BY  from  vector  Y,  calculated  at  step  5ii)  above. 

5.2.4  Numerical  Results 

In  Tables  5«1  - 5*4  we  show  results  of  Algorithm 

5.4  run  on  the  following  2 problems. 

Problem  I > 

y + Xy  - 0 

y(0)  = y ( 1 ) = 0 

2 2 

Eigenvalues  are  = (n+1)  n , n=0,l,... 

Problem  II 1 

y ♦ (A  - x|x|)y  = 0 
y(-l)  = y(l)  *■  0 

a 2.4625B070 
X2  a 22.2077749 


I 
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For  each  problem,  TOL  was  set  at  10"^.  The  columns 
are  Interpreted  as  follows.  N Is  the  site  of  the  matrix 
problem,  |e|  * I A - M,  lee8tl  the  error  estimate  cal- 
culated from  the  expression  shown  above  in  Step  3 of 
Algorithm  5*4,  |ecl  * A * egst  is  the  corrected  eigen- 
value resulting  from  Algorithm  4.6. 


Table  5*1 

Results  for  Problem  I,  X 
With  Algorithm  5*4 


N 

Ie| 

'Sst1 

»•«* 

8 

9.839D-4 

9.592D-4 

2.471D-5 

16 

6.122D-5 

6.091D-5 

3.127D-7 

32 

3.825D-6 

3.817D-6 

4.573D-9 

64 

2.388D-7 

2.387D-7 

6.899D-11 

Table  5.2 

Results  for  Problem  I, 
With  Algorithm  5*4  1 


N 

Ie| 

Kst1 

■•o' 

8 

7.479D-1 

1.011D0 

2.633D-I 

16 

4.516D-2 

4.313D-2 

2.031D-3 

32 

2.794D-3 

2.763D-3 

3.166D-5 

64 

1.742D-4 

1.737D-3 

4.698D-7 

128 

1.088D-5 

1.087D-5 

7.266D-9 

256 

6.801D-7 

6.798D-7 

2.932D-10 

II 


Table  5.3 
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Results  for  Problem  II,  With  Algorithm  5.4 


N 

lei 

''.St1 

lecl 

8 

2.278D-4 

2.232D-4 

4.597D-6 

16 

1.432D-5 

1.383D-5 

4.85  D-7 

32 

9.07D-7 

8.629D-7 

*4.  D-8 

Table  5*4 

Results  for  Problem  II,  With  Algorithm  5.4 

N 

lei 

I'.st' 

,ec' 

8 

1.871D-1 

2.524D-1 

6.529D-2 

16 

1.130D-2 

1.079D-2 

5.064D-4 

32 

6.991D-4 

6 . 911D-4 

8.02  D-6 

64 

4.363D-5 

4.346D-5 

1.7  D-7 

128 

2.77  D-6 

2.720D-6 

*5.  D-8 

256 

2.2  D-7 

1.701D-7 

*5.  D-8 

* Correct  to  all  available  digits 

The  discrepancy  between  |e|  and  |eegtl,  in  Tables 
5.3  and  5*4,  is  attributed  to  the  fact  that  and  are 
known  only  to  9 significant  figures  for  Problem  II.  The 
tables  show  us  that  Algorithm  5.4  can  be  used  effectively 
to  solve  the  S-I  problem.  The  error  estimate  appears  quite 
good  and  serves  well  as  a stopping  criterion  for  the  algo- 
rithm in  these  cases.  We  also  note  that  the  final  result 
of  the  algorithm  will  be  much  better  than  TOL  required 
because  the  corrected  value,  l e c | . gives  us  several  more 
significant  figures  of  accuracy. 
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5.3  An  Adaptive.  Iterative  Method 

We  recall  from  Section  5*1  that  an  adaptive  method 
is  one  in  which  the  mesh  points  are  chosen  to  depend  on  the 
solution,  (y,M  • This  obviously  implies,  in  general,  that 
the  mesh  will  be  non-uniform,  Thus  the  need  for  the  finite  - 
difference  method  of  Section  4.3  becomes  apparent.  If  the 
solution  were  known  beforehand,  we  could  choose  the  mesh 
points  appropriately  and  then  use  Algorithm  4.5  to  solve 
for  an  approximate  solution  on  our  optimal  mesh.  This  is, 
of  course,  not  possible  and  so  an  iterative  strategy  must 
be  devised  which  alters  the  mesh  to  conform  to  a converg- 
ing solution. 

Denny  and  LandisLl4]  solve  the  problem 

u + P(u,y)u  + Q (u , y ) = 0 
u ( 0)  = 0,  u ( 1 ) ---  1 


with  finite-difference  techniques.  They  obtain  successive 
distributions,  [y!k^},  of  the  independent  variable  by 
making  the  truncation  error,  T)  , go  to  zero.  To  accomp- 
lish this  they  must  continually  compute  dl\/<iy^  ^ , 
dTi/dyi'  and  dTi/dyi+l*  Eigenvalue  problems  are  not 
considered. 

Gary  and  Helgasonl.19] , who  do  solve  eigenvalue 
problems  with  finite-difference  methods,  require  that  the 
user  supply  functions  s(x)  and  s (x),  which  describe  the 


non-uniform  mesh.  They  then  transform  the  resulting  matrix 
problem,  A + XB,  into  the  standard  A ♦ XI  and  solve  with 
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the  QR  algorithm. 

Returning  again  to  Lentini  and  Pereyra£26],  we  find 
the  following  strategy  for  solving  first  order  systems 
adaptively.  In  each  interval  (x^,x.^)  they  add  J uniform- 
ly distributed  points,  where 

j - (truncation  error  in  the  interval) 

(that  interval's  alloted  error) 

They  note  that  unless  the  alloted  error  is  chosen  judi- 
ciously and  altered  throughout  the  process,  too  many  points 
will  be  introduced  into  the  mesh  early.  This,  of  course, 
defeats  the  purpose  of  an  adaptive  method. 

It  is  our  belief  that  the  optimal  strategy  would  be 
to  add  only  one  point  to  the  mesh  at  each  iteration.  A 
technique  like  the  one  described  by  Malcolm  and  SimpsonL2fi] 
for  numerical  integration  would  then  be  needed.  At  each 
iteration,  a point  would  be  added  to  the  interval  in  which 
the  truncation  error  is  largest.  An  elaborate  system 
of  sorting  or  a heap  type  data  structure  would  be  neces- 
sary, and  the  method  would  undoubtedly  become  very  inef- 
ficient. 

We  present  a compromise  to  these  various  techniques. 
At  most  one  point  is  introduced  per  interval,  specifically 
at  the  midpoint,  but  all  intervals  not  meeting  a certain 
criterion  are  refined.  Points  are  not  removed  from  the 
mesh  once  they  enter.  It  is  believed,  by  Lentini,  Pereyra, 
and  others  that  removal  of  points  can  prove  to  be  quite  un- 
stable and  in  addition  can  produce  very  coarse  meshes. 


n 
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*5 . 3 - 1 Refinement  Strategy 

Our  strategy  for  refinement  is  described  as  follows. 

Let  x.  < x.+  1 < ...  < x.  . be  a sequence  of  consecutive 
J J*  X Jt  K 

points  in  the  current  mesh  which  do  not  meet  some  given 
criterion.  (Such  criteria  will  be  discussed  below.) 

The  mesh  currently  looks  likes 


n 


H 1 

xj  xj+i 


* 1 1 

xj«k-l  xj*k 


Figure  5.1 

Mesh  Before  Refinement 


Points  [x.  ]P^+1  are  now  added  to  the  mesh  producing  the 

X 1 — J 

following  configuration. 


X j+kM 

-H  t i - r 


X . . . , X . 

Jl+k-l  j+k 


Xj  x.i+i 

t i -T  i-  4-4- 

xj  xjn 


Figure  5*2 

Mesh  After  Refinement 


Each  x^  is  at  the  midpoint  of  the  interval 

This  is  done  for  each  interval  of  points  not  satisfying 

the  criterion.  The  only  exception  is  if  x ^ or  Xj+k  is  an 

endpoint,  i.e.  x.  - a or  x-  . = b.  In  this  case  no  point 
J J * 

is  added  outside  the  interval  l_a,b]. 


DEFINITION  s When  a point,  x.,  is  refined  by  adding 

J 

x.  and  xj^^  as  described  above,  we  say  it  has  been  split. 


rt 


An  important  consequence  of  our  refinement  strategy 
is  the  effect  that  the  splitting  of  a point  has  on  the  vec- 
tor t(v) • We  see  from  the  following  theorem  that  the  de- 
sired result  does  in  fact  take  place. 


THEOREM  5 . 5 » When  x.  is  split  by  the  above  strategy 

J 

and  replaced  by  the  points  x^,  xj2*  an<*  xjv  then 


IT(x.)|  > £ IT(x..)| 

.]  k=l  JK 


(^.3.1) 


provided  h = max(h.)  is  sufficiently  small  and  “C(x.)  / 0. 

i 1 J 

Before  proving  the  theorem  we  make  the  following 

observations.  The  restriction  on  T(x.)  not  being  zero  is 

J 

certainly  reasonable,  for  no  strategy  should  attempt  to 
split  such  a point.  The  restriction  we  place  on  h,  that  of 
being  sufficiently  small,  is  not  a real  restriction  in  most 
practical  problems.  It  is  necessary  to  insure  that  the 
higher  order  derivatives  at  the  two  new  points,  x.^  and 
x^,  are  not  much  larger  than  those  at  x . . If  such  were 
the  case  then  h.  and  h.+,  would  need  to  be  smaller  before 

u ^ 

the  theorem  is  valid. 


Proof  of  Theorem  Since  we  are  assuming  all 
along  that  the  solution  y has  as  many  continuous  deriv- 
atives as  we  need,  we  know  that  for  all  € > 0,  h can  be 
made  small  enough  to  ensure  that 


(m)  - W 

|ylm,(xji)  - y ( x j ) | < € 

|y(m^(x.^)  -y(m^(Xj)|  < €,  m=5i6, . . . 

We  now  look  at  the  physical  changes  brought  on  by 
the  splitting  of  x.. 


Xj-1 


h.i+l 


Figure  5»3 


Xj+1 


Point  x.  Before  Splitting 

J 


h.  h . 


M jl 


Xj+1 


Figure  5.4 


Point  x.  After  Splitting 

We  see  from  Figures  *5 . 3 and  5.4  that  at  the  point  x^, 
or  x , the  only  change  is  that  both  h.  and  h are 
halved.  From  equation  (5.2.3)  for  the  truncation  error, 


we  can  now  say  that 


I t(x.2)|  < IX  (Xj)|  / 32 


and  in  fact  if  h^  - h^^  we  have 


(5-3.2) 


I t(x.2)|  < I T (Xj)  | / 64. 


(5.3*3) 


The  two  new  points,  x.^  and  x.^.  are  located  at  the  mid- 
points of  existing  intervals  and  thus  contain  only  even 
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terms  in  their  truncation  error  expansions. 

If,  in  the  original  configuration  (Figure  5*3) » we 
had  had  h . = h .+, , then  the  theorem  is  immediate  because 
the  points  and  x^  will  also  satisfy  inequalities  like 
(6.3-3).  i.e. 

I X (XjjH  < I t(x.)|  / 64 
| t (xj3)|  < | t(x.)|  / 64 

provided,  of  course,  that  the  derivatives  at  x^  and  x^ 
are  not  too  large . Thus 

|t(xjR)l  < 31  T (Xj)  |/64  < IT(x.|. 

If  the  original  configuration  had  h^  / h^+^  then 
we  know  that  the  truncation  error  at  x.  contains  odd  terms 

J 

and  in  particular  an  h**  term.  The  two  new  points  contain 
h^  terms  as  their  largest.  They  also  have  h.'s  which  are 

•J 

half  of  the  original.  Thus  coupled  with  (5-3*2)  it  is 
obvious  that  (6*3.1)  holds  again.  (Again  provided  that  the 
derivatives  are  sufficiently  smooth.)  Thus  the  theorem 
is  proved. 


In  addition  to  the  smaller  truncation  errors  at 
these  points  we  should  note  that  at  xj_^  and  xj+i  one 
the  h's  is  halved  and  thus  IT(x^  1)|  and  IT(Xj+1)| 
are  also  decreased  by  splitting. 


While  the  actual  amount  of  decrease  seems  quite 
difficult  to  compute,  it  is  obvious  from  expression  (4.3.21) 
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that  the  act  of  splitting  will  make  smaller.  We  also 

4 

know  from  Theorem  4.19  that  0(h  ) and  thus  several 

acts  of  splitting,  which  decrease  h,  are  guaranteed  to 
decrease  I . 

5.3.2  Refinement  Criteria 

Now  that  the  mechanics  of  splitting  a point  have 
been  described  we  must  have  some  criterion  for  deciding 
which  points  should  be  split.  Several  strategies  have 
been  tried  and  we  present  each  here  along  with  some  numer- 
ical results  comparing  them.  Based  on  these  results  and  on 
the  analysis  of  each  strategy  a "best"  strategy  for  solving 
our  problem  is  chosen.  By  "best"  we  mean  here  the  strategy 
which  most  efficiently  (fewest  number  of  points)  computes 
the  eigenvalue  to  the  desired  accuracy.  It  is  probable 
that  if  something  other  than  the  eigenvalue,  e.g.  the 
eigenfunction  only  or  the  expection  of  operators,  were  the 
primary  objective,  another  strategy  may  be  preferred. 

5 . 3 . 2 . 1 Equidistribution  Strategy 

Pereyra  and  Sewell[33J  suggest  a strategy  in  which 
X(x^),  i-1.2,...,N,  is  approximately  constant.  This  is 
the  strategy  used  by  Lentini  and  Pereyra  to  solve  their 
first  order  systems  of  boundary  value  problems.  Imple- 
mentation of  such  a strategy  is  quite  easy.  We  simply  split 
any  point  x^  at  which  |T(x^)|  is  greater  than  some  toler- 
ance, EPS.  We  then  resolve  the  matrix  problem  on  the  new 
mesh  and  repeat.  We  know  from  Theorem  5.5  that  eventually 
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all  IT  (x^)|  will  be  less  than  EPS  . We  have  already  ar- 
gued that  this  process  will  also  decrease  |aA | . We  show 
some  numerical  results  using  this  strategy  after  all  strat- 
egies have  been  considered. 

5 . 3 . 2 . 2 Solution -Weigh ted  Strategy 

Expression  (4.3-21)  gave  us  an  error  estimate  for 
aA  . We  can  also  get  a less  accurate  estimate  by  the  fol- 
lowing reasoning. 

We  assume  now  that  an  accurate  eigenvalue  is  our 
objective  and  that  the  error  in  the  eigenfunction  is  un- 
important compared  to  the  error  in  the  eigenvalue.  Re- 
calling now  that 

Ay  - XBy  = T(y)  (5-3-4) 

we  let  A = A ♦ a X and,  by  our  assumption,  we  let  y - Y. 
Substituting  into  (5*3.4)  and  multiplying  by  Y*  we  get 

Yt(AY  - ABY)  - AA Y^BY  » Y*  T(y)  . 

The  first  term  is  zero  by  expression  (4.3.2)  and  so  we 
see  that 

aA  » - Y 1 t(y)  / Y^BY . 

The  numerator  can  be  written  as 
♦ _ N+l 

- Yt  T(y)  -•  - £ Y.  T(x. ) , 
i-0 

and  so  from  this  estimate  of  aA  it  is  obvious  that,  while 

n 

n 
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making  IT(x^)|  small  will  decrease  | , a better  strat- 
egy should  be  to  make  |Y^T(x^)|  small. 

So  our  second  strategy  is  to  split  any  point  for 

which 

|Yi  T(x^) | > EPS 

thus  leading  to  final  mesh  on  which 

I Y i X (xi> | a | Y j T (x j) I for  all  i,j. 

5 . 3 . 2 . 3 D"1  Solution-Weighted  Strategy 

Again  considering  the  error  estimate  (4.3.22) 

AX  « -Y^D”1D“1  t(y)  / YtD‘1D"1BY  (5.3-5) 

we  introduce  our  third  strategy.  The  numerator  of  (5.3*5) 
is 

N+l  , 

- L d7f  Y.  t (x.)  (5.3.6) 

i=0  11  1 1 

where  the  (d.^}  are  elements  of  D.  If  this  expression  can 
be  made  small  then  |^X | will  also  be  small. 

The  first  two  strategies  do  in  fact  make  (5.3*5) 
smaller.  However,  expression  (5.3*6)  implies  that  what 
should  actually  be  equidistributed  is  d 7?  Y^t(x^).  Thus 
our  third  strategy  is  to  split  any  point  x^  for  which 
|d7?  Y^  T(x^)|  > EPS.  This  will  lead  to  a mesh  on  which 

ldii  Yi  ^ (xi^  58  Yj  ^ for  a11  ifj* 

Implementation  of  this  strategy  is  slightly  more 

difficult  because  the  matrix  D”1  would  not  normally  be 

il 

0 

[ l— 


w 


„ 


L 


^ ^ 


105 

computed.  We  must  recall  from  Section  4.3  that  the  coef- 
ficients <Xq£»  a2i’  and  **2i’  tha*  make  up  d^,  must 

be  recomputed  for  each  new  mesh. 

In  Table  5-5  we  compare  the  3 strategies  by  using 
each  to  solve  several  problems.  The  value  of  EPS  is  held 
constant.  We  show  the  number  of  points  (N)  required  and 
the  error  in  the  eigenvalue  (|e|)  resulting  in  each  case. 

The  results  in  Table  5*5  are  but  a few  of  those  we 
computed  but  they  are  quite  typical  of  the  observed  phe- 
nomena. The  equidistribution  strategy  uses  the  most 
points  in  almost  every  case,  however,  the  errors  are  quite 
similar.  The  other  2 strategies  tend  to  be  indistinguish- 
able when  solving  the  regular  problems  of  this  chapter. 

This  can  best  be  explained  by  considering  an  example . 
We  look  at  the  problem 

y = - *-y,  y(o)  = yd)  = o 

which  generated  the  first  row  of  Table  5*5.  The  final 

mesh  of  78  intervals  consists  of  only  2 different  sizes 

of  h.  Either  h = 1/64  or  h - 1/128.  This  relatively 

smooth  mesh  produces  a matrix  (D  ) which  has  a 1 at  37 

of  its  79  diagonal  elements.  The  other  elements  are  either 

.2 

1.5  or  2.  Thus  the  inclusion  of  the  terms  d^  into  the 
splitting  strategy  has  very  little  effect  on  the  final  mesh 
when  the  eigenfunctions  are  well-behaved  like  those  of  this 
chapter. 


Table  5.5 

Comparison  of  Refinement  Strategies 
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Equidist 

SW 

D'1 

- SW 

q(x) 

N 

lei 

N 

Ie| 

N 

Ie| 

0 

158 

1.56D-? 

78 

6.43D-8 

78 

6.43D-8 

6cos( 2x) 

126 

9.23D-8 

102 

1.20D-6 

98 

1.98D-7 

l6c os ( 2x) 

66 

6.82D-6 

54 

2.84D-5 

58 

5.06D-6 

x|x| 

38 

5.86D-7 

28 

6.91D-6 

34 

1.19D-6 

To  make  a decision  on  which  strategy  to  use  we  look 
at  the  error  estimates  which  motivated  the  2 strategies. 
Using  the  same  problem  as  an  example , we  find  that  for  both 
strategies  |e|  » 6.430x10'®,  but  leestl  for  the  SW  strat- 
egy is  about  1.163x10“^  while  for  the  D'1  - SW,  I ee s ^ I * 

•8 

6.428x10'  . Similar  differences  are  observed  for  all  the 
problems  looked  at.  Thus  we  conclude  that  when  accurate 
computation  of  the  eigenvalue  is  our  primary  objective,  the 
D“*  solution-weighted  strategy  should  be  used.  This  is  the 
strategy  used  to  obtain  all  further  results  of  this  chapter 
and  of  Chapter  7. 

5.3.3  The  Adaptive.  Iterative  Algorithm 

We  now  present  the  adaptive,  iterative  algorithm. 

ALGORITHM  5.6 1 

1.  Set  up  a basic  uniform  mesh,  nQ,  on  Qa.b]. 

2.  Solve  the  matrix  problem  (4.3.2)  using  Algo- 
rithm 4.5. 


r 


r 


r * 


ii 
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3.  Compute  error  estimate  (4.3*22)  and  make  one 
correction  using  Algor tihm  4.21. 

4.  If  |5X|  < TOL,  STOP. 

5.  Refine  mesh: 

i)  Split  all  points  x^  for  which 

Id if  Y.  T(Xi)|  > EPS. 

(See  below  for  a discussion  of  how  to 
compute  EPS  from  a given  TOL.) 
ii)  Interpolate  values  of  Y^  at  new  mesh 
points  to  get  starting  right-hand  side 
for  next  iteration. 

6.  GO  TO  2. 

In  Step  5 we  mention  a relationship  between  TOL 
(the  accuracy  desired  in  A)  and  EPS  (the  equidistribution 
level  at  each  point  of  the  mesh) . It  is  unreasonable  to 
expect  the  user  of  this  algorithm  to  specify  EPS,  though 
he  must  specify  TOL.  The  EPS  required  by  the  splitting 
strategy  can  be  computed  for  each  iteration. 

We  want 

iam  . i*yvfosn  < T0L. 

| YtD-iD_iBY| 

This  will  be  accomplished  if 


|YtD'1D”1T(Y) | 

< TOL* |YtD_1D"1BY| 

(5.3.7) 

|YtD"1D"1T(Y) | 

N+l  , 

< Z |d:f  Y.  T(x. ) | . 
i=0  11  1 1 

but 
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Since  each  |d^  T(x^)|  is  weighted  equally, 

(5.3*7)  will  be  satisfied  if 

,dIi  Yi  T(xi)l  < JU+2T* TOL* |YtD”1D*1BY| . 

Thus  if  we  let 

EpS  = t^.T0L-|yV1D-1BY|  (5*3*8) 

the  tolerance  will  be  met.  In  practice,  however,  we  found 
that  this  choice  of  EPS  is  too  small.  Using  it  the  actual 
|aX|  turns  out,  in  many  cases,  to  be  much  less  than  TOL 
thus  causing  extra  work.  We  attribute  this  to  at  least 
the  following  phenomena. 

1.  The  terms  of  Y^D'^D"^  T(y)  are  both  positive 
and  negative.  Thus  there  is  a considerable  amount  of 
cancellation  occurring  in  the  actual  aX . 

We  have  found  that  the  following  heuristic  gives 
better  results. 

a.  Count  the  number  of  times  POS  = Y^T(x^)  > 0 
and  NEG  = Y^T(x^)  < 0. 

b.  Instead  of  N+2  in  (5*3*8)  use  I POS  - NEG| . 

If  IPOS  - NEG | = 0 we  use  1. 

Our  reasoning  is  that  if  all  |d7?  Y^  *C(x^)|  = EPS  then  POS 
and  NEG  would  cancel  out  in  Y^D'^D-1  T(y)  and  only  the 
excess  one  way  or  the  other  would  affect  AX. 

2.  If,  on  iteration  j,  |d7?  Y^  T(x^)|  is  only 
slightly  larger  than  EPS,  we  find  that  splitting  nodes 
makes  it  much  less  on  iteration  ,1+1.  This  tends  to  make 


r 
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the  final  too  small. 

This  is  a difficult  occurrence  to  prevent  because 
there  is  no  way  to  foresee  the  exact  effect  splitting  will 
have  since  it  depends  on  all  points  of  the  mesh  and  their 
interaction. 

One  possible  solution  would  be  to  let 

EPS  ’ c'|p6S  - WEST  •TOi-irtD-1n-1BV| 

where  C is  some  constant  > 1.  Then  if  we  reach  a point  in 
the  iteration  where  |^A|  > TOL  but  all  IdT?  T(x^)|  < 

EPS,  we  decrease  G.  This  process  could  be  quite  ineffi- 
cient if  G were  not  decreased  b-  the  proper  amount. 

We  do  not  use  such  a "sliding"  EPS  in  our  test 
programs.  Rather  we  are  satisfied  with  an  actual  |aA| 
smaller  than  the  requested  TOL. 

5 . 3 . 4 Numerical  Results 

We  wish  to  compare  the  performance  of  Algorithm  5.^ 
with  that  of  the  uniform  mesh  method  of  Section  4.2.  'We 
tested  our  algorithm  on  each  of  the  following  problems 
with  the  results  shown  below.  After  the  algorithm  found 
its  optimal  N we  ran  the  Numerov  method  using  the  same  N 
and  the  2 closest  powers  of  2.  This  was  done  because  the 
non-adaptive  method  of  Section  5*2  is  restricted  to  N 


I 


being  a power  of  2. 

Problems  I and  II  are  described  in  Section  5*2.4. 
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Problem  III>  (Weber's  Equation) 

y (*  - x2)y  = 0 
y(0)  = y ( 1)  = 0 

XQ  a 10.1511640305 
as  247.0715002280 

Problem  IV t (Mathieu's  Equation) 

H 

y ♦ (X  - 2qcos2x)y  ~ 0 
y (0)  = y ( rr)  = 0 

n - 1 

XQ  a -.1102488 
X^  a 25.0208408 

q-.=  J 

Xx  a 3.2769220 
X3  a 16.2727010 
q = 8 

XQ  a -10.6053681 
X4  a 26.2209995 

While  the  results  shown  in  Table  5.6  do  display  some 
advantages  of  the  adaptive  method,  we  find  it  necessary 
to  point  out  that  the  method  does  not  perform  as  well  on 
regular  problems  as  on  the  singular  problems  of  Chapter  7. 

An  adaptive  strategy  seems  much  better  suited  to 
problems  in  which  the  solution  y(x)  has  singularities  or 
long,  exponentially  decaying  "tails"  rather  than  ones  in 
which  the  solution  is  more  well-behaved. 


Table  5.6 

Results  of  Algorithm  5.6 
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Problem 

Method 

N 

le| 

'eest' 

i.  k0 

adaptive 

78 

6.430D-8 

6.426D-8 

3.706D-11 

Numerov 

64 

2.391D-7 

2.387D-7 

7.006D-11 

78 

1.082D-7 

1.082D-7 

2.146D-11 

128 

1.492D-8 

1.492D-8 

1 . 258D-12 

I. 

adaptive 

136 

4.404D-5 

4.397D-5 

6.615D-8 

Numerov 

128 

6.115D-5 

6.108D-5 

7.109D-8 

136 

4.798D-5 

4.793D-5 

4.957D-8 

256 

3.821D-6 

3.819D-6 

1.130D-9 

u.  x0 

adaptive 

74 

*6.5  D-8 

2.225D-8 

*4.3  D-8 

Numerov 

64 

*6.7  D-8 

5.388D-8 

*1.4  D-8 

74 

*4.3  D-8 

3.015D-8 

*1.3  D-8 

128 

*1.5  D-8 

3.367D-9 

*1.2  D-8 

II. 

adaptive 

140 

9.368D-6 

9.306D-6 

*6.13  D-8 

Numerov 

128 

1.533D-5 

1.527D-5 

*6.18  D-8 

140 

1.073D-5 

1.067D-5 

*5.45  D-8 

256 

9.996D-7 

9.551D-7 

*4.45  D-8 

HI.  *0 

adaptive 

37 

1.719D-6 

1.714D-6 

*5-3  D-9 

Numerov 

32 

4.137D-6 

4.131D-6 

*6.3  D-9 

37 

2.314D-6 

2.312D-6 

*2.7  D-9 

64 

2.585D-7 

2.584D-7 

*1.  D-10 

III.  *4 

adaptive 

147 

4.007D-4 

4.007D-4 

1.18  D-8 

Numerov 

128 

2.333D-4 

2.329D-4 

4.197D-? 

147 

1.341D-4 

1.339D-4 

1.848D-7 

256 

1.458D-5 

1.457D-5 

*6.7  D-9 

Table  5.6  (continued) 


112 


Problem 

Method 

N 

lei 

'Set' 

le 

o' 

IV,  XQ 

adaptive 

70 

3.719D-7 

3.546D-7 

*1.7 

D-8 

q ■ i 

Numerov 

64 

2.601D-? 

2.426D-7 

*1.7 

D-8 

70 

1.868D-7 

1.696D-7 

*1.7 

D-8 

128 

*3.2  D-8 

1.518D-8 

*1.7 

D-8 

iv,  x. 

adaptive 

152 

1.03  D-5 

1.029D-5 

*1. 

D-8 

q - 1 

Numerov 

128 

2.38  D-5 

2.378D-5 

*2. 

D-8 

152 

1.20  D-5 

1.196D-5 

*1. 

D-8 

256 

1.4?  D-6 

1.487D-6 

*1. 

D-8 

IV,  x 

adaptive 

112 

*3.2  D-8 

6.888D-8 

*3. 

D-8 

q * 3 

Numerov 

64 

7.47  D-6 

7.414D-6 

*6.1 

D-8 

112 

8.23  D-7 

7.918D-7 

*3.3 

D-8 

128 

4.95  D-7 

4.642D-7 

*3.1 

D-8 

IV,  X 

adaptive 

128 

8.640D-6 

8.816D-6 

1.8 

D-7 

q - 3 

Numerov 

128 

7.100D-6 

7.283D-6 

1.8 

D-7 

iv,  xQ 

adaptive 

114 

5.5  D-7 

2.5HD-8 

*5. 

D-8 

q = 8 

Numerov 

64 

6.9  D-6 

6.827D-6 

*8. 

D-9 

114 

7.2  D-7 

6.8OOD-7 

*4. 

D-8 

128 

4.6  D-7 

4. 280D-7 

*4. 

D-8 

IV.  x4 

adaptive 

164 

2.057D-5 

2.050D-5 

*7. 

D-8 

q ■ 8 

Numerov 

128 

3.560D-5 

3.547D-5 

1.4 

D-7 

164 

1.322D-5 

1.317D-5 

*5. 

D-8 

256 

2.25  D-6 

2. 220D-6 

*3. 

D-8 

* Correct  to  all  available  digits 
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Even  with  these  concessions  we  find  that  in  half  of 
the  data  sets,  |e|  is  smaller  with  the  adaptive  method  than 
with  the  Numerov  method  using  the  same  N.  More  signifi- 
cantly we  find  that  if  a non-adaptive  method  like  that  of 
Section  5.2  were  used,  we  would  need  to  go  to  the  next 
highest  power  of  2 in  almost  every  case  to  get  equivalent 
or  better  accuracy. 

Another  result  of  Table  5.6  is  the  confirmation  of 
our  analysis  concerning  the  asymptotic  order  of  convergence. 
Theorem  4.19  assures  us  that  A - A * O(h^)  and  the  results 
of  Section  4.3.7  show  that  A + - A - O(h^). 

During  the  running  of  Algorithm  5*6,  N proceeds  from 
a value  of  8,  through  several  intermediate  values,  to  its 
final  value  as  shown  in  Table  5.6. 

If  |e|  and  |ecl  are  computed  when  N - 8 (call  them 
|e|p  and  leclg)  then  the  orders  of  convergence  can  be 
computed  by  solving  the  following  equations* 


I e | q lei 

( 1/8 )a  (I/N)a 


and  I ec I 8 

(1/8)* 


( I/N )* 


If  the  theoretical  results  mentioned  above  are  cor- 
rect then  a should  be  4 and  d should  be  6.  Using  Problem 
I as  our  example,  we  find  that  for  AQ  we  must  solve 
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and 


9.839xl0~4  _ 6.430xl0~8 

• (rr/8)a  (n/?8)a 


5.125xlO~5  = 3 ■ 706xl0~11 

(tt/8)4  (n/7fi)tf 


to  get  a = 4.23  and  3 = 6.21.  For  a " 4.05  and 
3 - 5.89*  Similar  values  are  obtained  from  the  other 
problems  when  enough  decimal  places  are  available  in  the 
exact  A. . 

Figures  5-5  and  5.6  show  the  distriDution  of  the 
selected  points  for  and  of  Problem  III.  We  see 

from  these  graphs  that  only  a few  different  sizes  of  h are 
used  to  solve  a regular  problem.  As  we  mentioned  above, 
however,  it  is  the  placement  of  the  points  which  makes  the 
adaptive  method  more  advantageous.  More  dramatic  step 
size  changes  are  evident  in  the  singular  problems  of 
Chapter  7. 

In  the  figures, the  points  along  the  x-axis  are  the 
IxjJ  placed  by  the  algorithm.  The  graph  itself  shows  the 
eigenfunction  computed  by  the  method.  In  Figure  5.6  only 
half  of  the  eigenfunction  is  shown.  The  function  is,  how- 
ever, symmetric  about  x - .5. 


I 


* 


Mesh  Distribution  for  Problem  III 


Chapter  6 


STARTING  VALUE  AND  RECYCLING  PROBLEMS 

6.1  Starting  Value  Problem 

In  the  algorithms  of  Chapters  4 and  5 it  was  assum- 
ed that  a good  estimate,  *°  the  des*red  eigenvalue, 

(or  A^),  was  always  available. 

In  some  applications  it  is  possible  to  get  such  an 
estimate.  It  may  be  obtained  through  theoretical  means 
by  comparing  the  given  problem  with  a known  problem,  or 
when  solving  a specific  problem  which  represents  a physi- 
cal process,  an  estimate  may  be  obtained  through  experi- 
mental results.  When  a good  approximation  is  supplied  it 
is  evident,  that  convergence  of  the  adaptive  algorithm  will 
occur  quite  quickly  and  the  computer  program  will  be  highly 
efficient. 

In  many  other  applications  the  user  cannot  supply 
an  estimate  to  the  eigenvalue.  Instead  he  knows  only  that 
he  wants  the  kth  smallest  eigenvalue,  A^,  or  equivalently 
the  eigenvalue  corresponding  to  yk(x),  the  eigenfunction 
with  exactly  k zeros  in  (a,b).  Our  objective  then  in  this 
section  is  to  develop  an  algorithm  for  computing 

6.1.1  Objectives  of  the  Algorithm 

We  believe  the  following  objectives  to  be  the  most 
important  in  developing  a starting  value  algorithm.  Thus 
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we  consider  each  in  our  discussion. 

1.  The  algorithm  must  be  efficient.  The  entire 
advantage  of  our  adaptive  algorithm  can  be  lost  if  there 
is  no  efficient  method  for  finding  a starting  value. 

We  argue  in  Chapter  4 that  our  non-uniform  finite- 
difference  method  gives  cubic  convergence,  while  the  shoot- 
ing method  algorithms  like  the  one  by  Bailey.  Gordon  and 
Shampine  (see  Chapter  3 for  a discussion)  give  only  quad- 
ratic convergence.  However.  Bailey,  et  al.,  claim  that 
their  algorithm  will  converge  to  the  proper  eigenvalue 
regardless  of  the  starting  value.  This,  of  course,  is  be- 
cause the  number  k (defined  above)  is  a parameter  in  the 
boundary  conditions.  Since  our  method  uses  inverse  iter- 
ation it  will  converge  to  the  eigenvalue  closest  to  the 
initial  value.  We  must  therefore  use  a procedure  which 
computes  a starting  value  without  substantially  reducing 
the  convergence  rate. 

2.  The  algorithm  should  be  no  more  accurate  than 
necessary.  This  idea  is  discussed  in  greater  detail  in 
Section  6.1.3*  We  show  there  that  the  starting  value  algo- 
rithm converges  much  more  slowly  than  the  adaptive  algo- 
rithm and  thus  we  do  not  want  it  to  operate  longer  than 


necessary. 


3.  The  algorithm  must  use  existing  techniques. 


By 


r 


this  we  mean  that  the  algorithm  should  be  designed  to  use 
as  many  of  the  same  procedures  as  the  main  adaptive  algo- 
rithm as  possible. 
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While  it  might  be  possible  to  use  a shooting  meth- 
od, for  example,  to  get  a starting  value,  it  would  mean 
introducing  a whole  new  set  of  procedures.  We  instead  find 
a starting  value  directly  from  the  matrix  problem  (4.3,2). 

6.1.2  Sturm  Sequences 

Let  us  assume  that  the  user  of  our  adaptive  S-L  al- 
gorithm supplies  us  with  the  index  of  the  eigenvalue  de- 
sired, i.e.  k is  the  input  when  is  being  requested.  We 
wish  to  show  that  an  appropriate  choice  for  an  approxima- 
tion to  kk  is  an  estimate,  A k 0»  "the  kth  matrix 

eigenvalue. 

We  recall  from  Chapter  4 that  the  eigenvalues  of  the 
S-L  problem  are  approximated  by  the  eigenvalues  of  the 
matrix  problem 

(A  - A B)Y  = 0.  (6.1.1) 

We  know  from  Theorem  4.19  that  for  every  eigenvalue,  Xk»  of 
the  differential  equation  there  exists  an  eigenvalue,  A, 
of  the  matrix  problem  (6.1.1)  such  that 

(A  - *k)  = 0(h4)  as  h - 0. 

To  get  convergence  in  the  other  direction  we  use 
Theorem  3.2  and  Lemma  4.1  of  129].  Because  our  discret- 
ization operator,  (4.3.1),  call  it  L^,  is  such  that 

lim  L.  u = L u 


(6.1.2) 


1 
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where  L is  the  normal  form  of  the  regular  S-L  operator,  we 
can  say  that  spectrum(Lh)  -*  spectrum(L)  as  h - 0. 

To  solve  singular  S-L  problems  we  truncate  the  in- 
terval (a,b)  to  a finite  interval  in  which  q is  continuous. 
In  order  for  the  limit  (6.1.2)  to  hold  as  h - 0 and  as  the 
truncated  interval  - (a,b)  we  must  require  that  q € C(a,b). 

We  recall  from  Chapter  2 that  singular  problems  can 
have  certain  regions  of  continuous  spectrum.  We  must  fur- 
ther restrict  our  class  of  problems  to  those  in  which  any 
continuous  spectrum  is  greater  than  the  point  spectrum. 

Using  these  facts  we  can  now  prove  the  following. 

THEOREM  6.1 i Let  k be  fixed  and  less  than  or  equal 
to  the  number  of  eigenvalues  in  the  point  spectrum  of 

(4.1.1)  and  let  any  continuous  spectrum  of  L be  greater 

than  the  point  spectrum.  Then  as  h - 0,  A;  of  the  system 

J 

(6.1.1)  converges  to  k.  of  (4.1,1)  for  all  j * 0,1,. ...k. 

J 

Proof i We  require  k to  be  fixed  because  we  do  not 
want  it  to  go  to  » as  h - 0.  The  theorem  holds,  however, 
for  any  choice  of  k in  the  proper  range. 

First  we  recall  from  Chapter  2 that  the  of  the 
point  spectrum  can  be  ordered  as  follows i 


• • • 
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Theorem  4.16  assures  us  that  for  h small  enough  the 
are  distinct.  In  addition  there  are  only  a finite  num- 
ber of  A^  and  thus  the  above  ordering  is  possible. 

We  now  prove  our  theorem  by  induction.  We  start 
by  showing  that  as  h - 0,  AQ  - XQ. 

We  know  that  there  exists  some  j for  which  A^  - XQ. 
Let  us  assume  that  Aj  is  not  the  smallest  matrix  eigen- 
value, i.e.  j / 0.  Then  AQ  - some  X^  > XQ. 

Define  X - Xn  = 6 > 0.  We  let  h be  small  enough 
m v 

so  that 

|Aj  - XQ|  < 6/2  and  IA0  - Xj  < 6/2. 

This  implies  that 

X + X 

XQ  - 6/2  < A j < XQ  + 6/2  = -2-g — - 

and  - X.  - 6/2  < A 0 < ♦ ‘/2 

which  in  turn  imply  that  Aj  < Aq.  This  is  a contradic- 
tion to  the  ordering  of  the  Aj  and  thus  Aq  - XQ. 

Let  us  now  assume  that  A:  * X.,  j=0,l, . . . ,k-l,  and 

J J 

show  that  Afc  - X^.  We  know  that  there  exists  a p such 


that  A *•  Xv.  Prom  our  induction  hypothesis  we  know  that 
p / 0,l,...,k-l.  We  assume  that  p / k,  i.e.  p > k.  It  is 
true,  however,  that  Av  **  some  < Thus  we  have 
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Letting  X - X = 6 > 0,  we  can  choose  h small 
m k 

enough  such  that 

IAp  - V < 4/2 

and  I Ak  - ^ml  < 6/2. 

Together  these  imply  that  A_.  < A,,  which  is  a contradic- 
tion  to  the  ordering  of  the  A:.  Thus  the  theorem  is 

J 

proved . 

This  theorem  tells  us  that  as  h becomes  small,  X, 

k 

will  be  approximated  by  the  kth  smallest  matrix  eigen- 

value. Thus  as  a starting  value,  Ak  Q.  for  Algorithm  5.6 
we  approximate  the  kth  eigenvalue  of  the  original  matrix 
problem. 

The  procedure  we  use  to  get  this  approximation  is 
based  on  the  Sturm  sequence  property  of  the  matrix  problem 
(6.1.1).  We  use  a bisection  method  described  by  Wilkinson 
[42]  to  locate  and  approximate  the  kth  eigenvalue. 

Define  the  following  sequence  of  functions  derived 
from  the  matrix  problem  (6.1.1). 

P0(u)  = 1 

p,(u)  - det  (A  - uB)  ^ 

= «n  - bu(u  - qj) 

P^(u)  --det  (A  - uB)^ 


n 
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• Ki  - ,,u<u  - ’i’1  pi-i(u) 

" la2.i-l  " 32,i-ll“  ' qi'^‘ 

ta0i  - 30i<u  - qi-l>>  pi-2(u) 
l"2| 3 p • • • pN • 


If  we  let  the  sequence  be  evaluated  for  some  uQ  and 
let  s(u0)  denote  the  number  of  sign  agreements  of  consec- 
utive members  of  the  sequence,  we  can  use  the  following 
result  proved  by  Wilkinson[42,  p.  300 J : 

The  value  s(uQ)  is  the  number  of  eigenvalues  of 

(A  - uB) Y = 0 which  are  strictly  greater  than  uQ. 

This  result  is,  of  course,  only  valid  if  the  se- 
quence pQ,  p^,  ...  , is  a Sturm  sequence.  However,  we 
know  from  the  proof  of  Theorem  4.16  that  there  does  exist 

an  h such  that  for  max  h.  = h,  the  sequence  is  a Sturm 

i 1 

sequence.  We  have  found,  however,  in  all  our  test  problems 
that  even  the  coarsest  mesh,  n * 8,  yields  real  eigenvalues 
which  are  reasonable  estimates  to  Therefore  in  most 

cases  there  does  not  appear  to  be  a need  to  decrease  h to 
ensure  that  we  have  a Sturm  sequence. 

One  check  must  be  made,  however,  which  may  lead  to 
increasing  N immediately.  If  the  user  asks  for  where 
k > N,  then  the  matrix  problem  cannot  supply  us  with  the 
proper  eigenvalue.  Thus  N must  be  increased  even  before 
the  adaptive  algorithm  runs.  This  check,  therefore,  is 
included  as  part  of  Algorithm  6.4. 
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6.1.3  Bisection  Algorithm 

We  use  the  following  algorithm  to  locate  the  kth 
eigenvalue  of  (6.1.1). 

ALGORITHM  6.2» 

1.  Let  bQ  > aQ  be  such  that  s(aQ)  > N-k  and 

s(bQ)  < N-k.  (We  then  know  that  Ak  lies  in 
the  interval  (a0,bQ).) 


2. 

Let  i --  1. 

3. 

Let  c i - (ai_1  ♦ bi_1)/2. 

4. 

Compute  p0(c^),  p^c^),..., 

pN(ci> • 

5. 

Compute  s(c^) . 

6. 

If  s(c^)  > N-k  take  a^  - c^ 

and  b^ 

If  s(c^)  < N-k  take  a^  = a^ 

^ and  b 

7. 

If  (b.  - a.  ) < f then  STOP 

with  A 

else  increment  i by  1 and  GO  TO  1. 


This  process  can  continue  until  Ak  0 is  as  accurate 
as  desired.  'We  see,  however,  that  the  algorithm  only 
yields  linear  convergence  to  A,,  because  (b.  - a.)  ■-* 

(b^  i ~ ai  j_)/2*  Therefore  we  must  only  continue  long 
enough  to  ensure  convergence  to  the  correct  eigenvalue. 

By  the  nature  of  inverse  iteration,  this  conver- 
gence is  only  made  certain  if  our  estimate,  A^  q»  is 
closer  to  A^  than  to  any  other  eigenvalue.  If  this  were 
not  the  case  we  would  get  convergence  to  the  closest,  and 
thus  the  wrong,  eigenvalue.  The  size  of  the  interval 
(a^,b^)  is  thus  dependent  on  how  close  the  eigenvalues  are. 
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We  estimate  A k by  c^^  = (a^  b^)/2.  We  know  that 

lies  in  the  interval  (a^,b^)  and  so 

I ^ k I *■  (b^  — a^)/2. 

If  we  require  also  that  b^  < Ak+1  and  a^  > Aj^  or  equiv- 
alently that  s(b^)  * N-k-1  and  s(a^)  3 N-k  then  we  see  that 

- Aj|  — “ a^)/2,  j X k 

and  thus  is  closest  to  Ak»  This  check  is  quite  easy 

to  implement  because  s(a^)  and  s(b^)  have  already  been 
computed  when  s(Cj^)  is  being  evaluated. 

One  final  point  to  be  considered  before  presenting 
our  algorithm  is  the  problem  of  finding  appropriate  aQ  and 
bQ.  WilkinsonC42]  notes  that  suitable  aQ  and  bQ  when  B = I 
are  + llAll^.  This,  of  course,  is  a consequence  of  the  well 
known  result 

N 

A)  < max  £ | a. . | 
l<i<N  j=l  1J 

(Varga[4l],  p.  17)  where  p(A)  is  the  spectral  radius  of  A. 

For  our  problem  aQ  and  bQ  can  be  computed  by  noting 
the  following  result. 

THEOREM  6.3»  Let  the  generalized  eigenvalue  problem, 
(A  - AB)Y  = 0,  be  such  that  B is  strictly  diagonally  dom- 
inant. Then 

N N 

max  | X.  | < {,  £ la.  A / ( I b.  . | - £ |b.,|)} 

k K j=l  1J  11  j=l  1J 

j/i 


r 


126 


where  i is  such  that  |Y^|  > |Y^|,  j / i. 


n 


1 .e . 


Proof i Let  k be  any  eigenvalue  of  the  generalized 
problem  and  let  Y be  an  eigenvector  corresponding  to  k. 
Let  Y.  be  the  largest  component,  in  absolute  value,  of  Y, 
|Y.|  > |Yj|  , j / i. 

Writing  the  ith  equation  we  have 

- Xbij)Yj  = 0 
N 

(b.  -A.  - a.  . )Y.  = >:  (a.  . - kb.  JY.. 

hi  i-i  i J iJJ 

U 


or 


So  by  the  triangle  inequality  we  find  that 


lbiiX  - aii>  lYl'  S laij  • Xblj'  |yj'- 


Dividing  by  |Y.J  yields 


lbiix  - *ii'  s Iai.i  - Xbio'  lyjl  / lyi 
- iA  laij  - Xbij'- 

Using  the  triangle  equality  on  both  sides  we  get 


or 


lbi.|  \k\ 


IMllb..| 


- Ia..|  < E |a.  J ► I X | E | b.  . | 


j/i 

j/i  J 


ij 

N 

E 

j-1 


j/i 


1J 


a 


ij1. 


Since  b.  ~ | b. . | - L I b,  . | > 0 
1 “ j/i  lJ 


we  have  finally 


and  the  theorem  is  proved. 


Since  we  do  not  know  which  component  of  Y is  largest 
before  Y is  found  we  use 

N 

t max  ( £ |a.  -I  / b. } (6.1. 3) 

i j*l  1J  1 

for  aQ  and  bQ . 

6.1.4  Starting  Value  Algorithm 

We  now  present  the  starting  value  algorithm. 

ALGORITHM  6.4t  let  N •-  8.  (arbitrary) 

1.  If  k > N then  N = 2N  and  repeat  1. 

2.  Set  up  matrices  A and  3. 

3.  Compute  aQ  and  bQ  from  (6.1.3). 

4.  Compute  A^  Q by  Algorithm  6.1  with  the  fol- 
lowing exception i 

At  Step  7 we  have 

7.  If  s(a^)  - N-k  and  s(b^)  --  N-k-1 
then  6 = (b^  - a^)/2  and  STOP  with 

Ak  o " ci+l  else  and  G0  T0  3* 

It  should  be  noted  that  when  Algorithm  6.4  is  com- 
pleted the  matrix  A - Ak  QH  is  exactly  the  one  we  require 
to  execute  Algorithm  6.6.  Thus  all  our  objectives  of 
Section  6.1.1  have  been  satisfied. 


f ' 

I Li 
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6.1.5  Numerical  Results 

We  present  here  2 examples  of  the  starting  values 
found  by  Algorithm  6.3.  The  problems  chosen  are  Problems 
I and  III  of  Chapter  5»  These  problems  were  selected  be- 
cause many  exact  eigenvalues  are  known  and  our  algorithm 
can  be  well  tested. 

In  Tables  6.1  and  6.2  we  show  several  eigenvalues, 
A^,  followed  by  the  estimate,  EV,  from  Algorithm  6.4.  The 
final  column  shows  the  number  of  bisections  needed  to  get 
the  value . 

In  each  case,  the  estimated  value  is  closer  to  the 
ith  eigenvalue  than  any  other.  We  also  note  that  as  i 
gets  larger,  the  number  of  iterations  decreases.  This  is, 
of  course,  due  to  the  fact  that  in  these  problems  the 
eigenvalues  are  more  greatly  separated  as  i gets  larger  and 
thus  not  as  many  bisections  are  necessary  to  ensure  con- 
vergence to  the  proper  eigenvalue. 


Table  6.1 

Results  of  Algorithm  6.3  on  Problem  I 


i 

*i 

EV 

# bisect 

0 

9.870 

6.0 

10 

1 

39.478 

42.0 

10 

2 

88.826 

84.0 

9 

3 

157.914 

168.0 

8 

4 

246 . 740 

264.0 

8 

5 

355.306 

360.0 

8 

TO 

1194.22 

1200.0 

7 

19 

3947.84 

3942.0 

6 

Table  6.2 

Results  of  Algorithm  6.3  on  Problem  III 


i 

X. 

X 

EV 

# bisect 

0 

10.151 

6.001 

10 

1 

39.799 

42.006 

10 

1 

2 

89.154 

84.012 

9 

3 

158.244 

168.024 

8 

4 

247.072 

264.038 

8 

5 

355.638 

360.051 

8 

10 

1194.555 

1200.172 

7 

1 

19 

3948.175 

3744.535 

6 
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6.2  Recycling  Problem 

We  now  consider  a problem  closely  related  to  the 
starting  value  problem  of  Section  6.1.  We  call  it  the 
Recycling  Problem  for  reasons  which  will  be  obvious. 

Recalling  Algorithm  5.6,  the  adaptive,  iterative 
algorithm,  we  note  that  after  each  iteration  of  the  algo- 
rithm we  must  set  up  a new  matrix  in  order  to  begin  the 
next  iteration.  This  requires  a value  of  A (call  it 
Ak  ^ for  the  ith  iteration  of  A^)  to  set  up  the  matrix 

A “ ^k,iB* 

A first  inclination  is  to  use  the  value  of  A found 
after  iteration  i-1  was  completed.  We  see  from  the  follow- 
ing numerical  example,  however,  that  this  can  lead  to 
incorrect  answers. 

The  problem  we  solve  is  a singular  S-L  problem: 

y"  , (-  J - My 

y ( 0 ) = y(M  * 0 

p 

having  eigenvalues  *n  - - l/(4n  ),  n~l,2,...  Such  a prob- 
lem is  much  more  likely  to  converge  to  a wrong  eigenvalue 
than  is  a regular  problem. 

Ve  attempt  to  compute  k - -.0625  using  the  tech- 
nique described  above.  Table  6.3  shows  the  values  of 
and  K as  N increases.  Recall  that  X is  the  corrected  eig- 
envalue. Algorithm  6.4  was  used  to  compute  a starting 
value  of  A -.03721.  Convergence  to  the  wrong  eigenvalue 
: r;  has  occurred. 

I 1 . 

i n 

u 
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Table  6.3 
Recycling  of  A^ 


N 

Ai 

X 

starting  value 

-.03721 

8 

-.03890 

-.0393^ 

18 

-.04101 

- . 04474 

25 

-.05910 

-.25541 

30 

-.1454? 

-.22907 

42 

-.23886 

-.29142 

51 

-.24418 

-.24553 

55 

-.24699 

- . 24704 

During  early  iterations  of  Algorithm  5*6,  when  the 
mesh  is  coarse,  the  kth  matrix  eigenvalue  may  be  quite  dif 
ferent  from  the  exact  It  may  even  be  closer  to  some 

other  X y This  is  seen  in  the  numerical  example  of  Table 
6.3.  This  phenomenon  is  especially  prevalent  when  higher 
eigenvalues  are  desired. 

A second  possible  solution  is  to  run  the  starting 
value  algorithm  before  each  iteration  of  Algorithm  5*6. 
This  assumes  no  knowledge  of  the  eigenvalue.  In  later  it- 
erations we  can  certainly  expect  that  the  recycled  value 
will  be  quite  close  to  the  desired  eigenvalue.  Thus  this 
second  strategy  would  tend  to  be  quite  inneficient. 

What  we  need  then  is  a procedure  to  check  A^  ^ ^ 
for  validity,  i.e.  whether  it  falls  within  ( A^ Aj^) 
and  whether  it  is  closest  to  A^.  If  these  conditions  are 
not  met,  the  procedure  must  produce  a value  for  which  they 
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are . 

We  claim  that  the  following  algorithm  will  fulfill 
these  requirements.  In  Steps  4.1  and  5.1  the  number  & is 
used.  The  reader  is  referred  back  to  Algorithm  6.4  to  find 
that  6 a (b^  - a^)/2,  i.e.  half  the  length  of  the  inter- 
val found  to  contain  Ak  Q.  The  algorithm  attempts  to 
show  that  the  interval  (Ak  Ak  + 6)  lies 

within  ( Ak_x.  Ak4l).  If  so  then  Ak>i-1  can  serve  as  the 
starting  value  for  iteration  i.  If  this  cannot  be  shown 
then  a new  value  must  be  found. 

ALGORITHM  6.5:  Let  Av  ; , be  the  matrix  eigenvalue 

K»  1-1 

computed  by  iteration  i-1  of  Algorithm  5-6. 

1.  Set  up  the  matrix  A - Ak.  ^ for  iteration  i. 

2.  Compute  sk  = s(  A<.i-i>- 

3.  If« 


= N-k-1 

GO 

TO 

4 

<A<An> 

= N-k 

GO 

TO 

6 

^-Ak,i-1 

€ 

< 

» 

H 

l 

& 

< 

< N-k-1 

GO 

TO 

6 

lAk.i-l 

€ 

<Akn.-)l 

> N-k 

GO 

TO 

7 

lAk,i-l 

( 

r—' 

1 

<r 

• 

S 

i 

4. 

4.1.  Let  bQ  = Ak>i_1. 

4.2.  Let  aQ  - bQ  - 6 . 

4.3*  Compute  s(aQ). 

4.4.  If« 

s(a0)  •-  N-k  then  STOP  and  recycle 

(a0  * V/2* 
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s(aQ)  N-k-1  then  bQ  - aQ  and  GO  TO  4.2. 

s(aQ)  > N-k  then  run  Algorithm  6.4. 

5- 

5*1*  Le  t aQ  Ak  ^ ^ * 

5.2.  Let  bQ  = aQ  + 6 . 

5.3.  Compute  s(bQ) . 

5.4.  Ift 

s(bQ)  - N-k-1  then  STOP  and  recycle 

U0  + b0V2. 

s(bQ)  - N-k  then  aQ  - bQ  and  GO  TO  5.2. 

s(bQ)  < N-k-1  then  run  Algorithm  6.4. 

6. 

6.1.  Let  bg  ~ ^^k  i_i • 

6.2.  Compute  aQ  from  (6.I.3). 

6.3.  Run  Algorithm  6.4. 

7. 

7.1.  Let  a0  - Akfl_r 

7.2.  Compute  bQ  from  (6.1. 3). 

7»3*  Run  Algorithm  6.4. 

In  the  early  iterations  of  Algorithm  6-6  we  may  find 
that  Algorithm  6.4  is  run  quite  often  at  Steps  4.4,  5.4, 

6.3,  or  7«3  of  Algorithm  6.5.  However,  in  later  iterations 
when  Ak  becomes  a good  approximation  to  Ak  we  find 
that  at  either  Step  4.4  or  5.4  the  algorithm  stops  the 
first  time  through.  This  signifies  that  . is  an  ade- 

quate  approximation  to  Ak< 


Using  this  algorithm,  we  reran  the  example  above. 

The  results  are  shown  in  Table  6.4.  The  column  headed 
"step  shows  the  path  taken  through  Algorithm  6.5  to 
obtain  each  value. 

This  example  shows  that  the  value  computed  at  iter- 
ation i can  indeed  be  a poor  estimate  of  the  desired  eigen- 
value. We  observe  how  Algorithm  6.5  reacts  to  update  that 
value  so  that  convergence  to  the  correct  value  is  assured. 
It  also  shows  that  in  higher  iterations  the  recycled  value 
can  be  an  adequate  estimate.  In  the  example  we  see  that 
from  N = 44  to  N = 54,  ^ ^ is  found  to  be  within  the 

proper  interval. 


Table  6.4 

Recycling  With  Algorithm  6.5 


i 

N 

A 

r 

step  # 

starting  value 

-.03721 

1 

8 

-.03890 

-.03934 

- 

2 

18 

-.1586? 

-.1328? 

6 

3 

25 

-.18042 

-31-987 

4 

4 

29 

-.25238 

-2.5246 

5 

5 

40 

-.06215 

-.0643? 

7 

6 

44 

-.06213 

-.06222 

5 

7 

48 

-.06213 

-.06214 

5 

8 

49 

-.06213 

-.06214 

4 

9 

51 

-.06231 

-.06232 

4 

10 

53 

-.06241 

-.06241 

4 

11 

54 

-.06241 

-.062 41 

4 

Chapter  7 
SINGULAR  PROBLEMS 

7.1  Introduction 

The  next  topic  we  discuss  is  that  of  singular  S-L 
problems.  Some  of  the  theory  of  the  singular  problem  is 
presented  in  Chapter  2.  In  the  present  chapter  we  show 
the  additional  capacity  of  our  adaptive  algorithm  to  solve 
the  two  most  common  classes  of  singular  problems.  ^ 

’We  recall  the  Sturm-Liouville  system  written  in 
Liouville  normal  form. 

y"  + (*•  - q (x)  )y  - 0 

y(a)  = 0 , y(b)  = 0 

For  problems  in  the  normal  form  we  find  that  singular 
problems  generally  fall  into  the  following  two  classes. 

a.  The  interval  (a,b)  is  infinite,  i.e.  a -•  -*>  or 
b - <*>  or  both.  The  boundary  conditions  of  such  a problem 
are  of  the  form  y(x)  * 0 as  x - », 

b.  The  function  q(x)  is  discontinuous  at  a or  b or 

both. 

Our  adaptive  method  is  equipped  to  handle  these  types  of 
problems . 


Both  types  of  problem  are  solved  by  truncating  the 
interval  (a,b)  and  solving  the  regular  S-L  problem  defined 
on  the  truncated  interval.  By  restricting  the  class  of 
problems  as  indicated  in  Section  6.1.2.  i.e.  q € C(a.b) 
and  any  continuous  spectrum  is  greater  than  the  point  spec- 
trum. we  are  assured  that  as  h -*  0 and  the  truncated  inter- 
val - (a.b),  the  eigenvalues  and  eigenfunctions  of  the 
truncated  problem  will  converge  to  those  of  the  singular 
problem. 

7.2  Infinite  Domain 

Existing  techniques  for  solving  the  singular  Sturm- 
Liouville  problem  with  an  infinite  domain  generally  lie  in 
two  categories.  First  is  the  technique  of  mapping  the  in- 
finite interval  onto  a finite  interval.  For  example,  from 
Chapter  3 we  recall  that  Bailey,  Gordon,  and  Shampine[l] 
transform  all  problems,  regular  or  singular,  to  the  inter- 
val (-1,1). 

The  second  technique  is  to  restrict  the  domain  to 
some  finite  interval,  i.e.  instead  of  (a,00),  we  solve  the 
problem  on  (a,R),  imposing  the  boundary  condition  y(R)  = 0, 
where  R is  some  appropriately  chosen  large  number.  Since 
we  have  already  assumed  that  all  problems  are  transformed 
to  normal  form  we  do  not  wish  to  require  even  further 
transformation.  Thus  we  choose  the  second  technique  as  the 
basis  for  our  method. 

The  question,  and  one  which  cannot  be  answered  with- 
out prior  knowledge  of  the  solution,  y(x),  is  how  large  to 
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... 


make  R.  If  R is  too  small  we  may  truncate  a significant 
portion  of  the  eigenfunction  thus  greatly  affecting  the 
eigenvalue.  If  R is  too  large  we  will  have  introduced  ex- 
tra mesh  points  which  have  no  effect  on  the  eigenvalue. 
This  would  be  counter-productive  to  the  goals  of  our  adap- 
tive method. 

To  be  truly  adaptive  our  algorithm  must  be  capable 
of  adjusting  R.  In  general  we  find  that  the  optimum  value 
of  R depends  not  only  on  the  problem  but  on  which  eigen- 
value is  desired  within  a given  problem.  Fischer  and 
Usmani[l?]  have  investigated  the  infinite  domain  with 
regard  to  two-point  boundary  value  problems.  We  use  some 
of  their  ideas  here. 

We  start  by  assuming  that  some  value  of  R is  avail- 
able. Some  knowledge  of  the  problem  may  make  a reasonable 
guess  at  R possible.  If  not,  the  algorithm  uses  some  pre- 
assigned value.  We  also  assume  that  b = * while  a is  fi- 
nite. All  results  are  easily  adapted  to  the  case  in  which 
a 3 - and  we  show  some  numerical  results  for  each  case 
in  Section  7.4. 

When  the  matrix  problem 

(A  - A B)Y  = 0 (7.2.1) 

i 8 solved  on  the  interval  (a,R),  we  make  the  assumption 
that  = 0,  however,  / 0*  Thus  the  actual  sol- 

ution (*,y)  of  the  differential  equation  satisfies 


i 


(A  - XB)y  = T(y)  - eR 
where  eR  = (0,0, . . . ,€R)*  with 

(R  * «2Ny<xIH-l)  * 82t/’(xN»l>  * 0 • 

This  is  seen  by  looking  at  the  last  equation  of  the  system 

(7.2.2), 

I0Ny(xN-l>  * 2y<V  * !,2Ny(xNU)  * S0Ny"(xN-l)  * 

SlNy  *XN*  * 82Ny  *xN+l'  = 

If  we  now  repeat  the  derivation  of  equation  (4.3*21) 
of  Section  4.3  starting,  however,  with  equation  (7.2.2) 
rather  than  with  (4.3.9)  we  find  that 

aX  = - YtD"1D“1(  t(y)  - iR}AtD-1D-1B7. 

If  we  assume  that  AXAy  is  negligible  we  get 

yW1  r(y)  - YtD“1D"1eR 

AX  * - r—. 5 — -.  (7.2.3) 

ytt»"1d*-lby 

Our  D”1  solution-weighted  strategy  is  based  on  the 
expression  (4.3.22).  We  found  that  we  should  require  that 

ldIi  T(xi)l  - ldi!  Yj  T<xj>'- 

We  now  see  from  (7*2.3)  that  the  site  of  should 

also  be  monitored. 

Since  eR  contains  zeros  everywhere  but  the  last  row. 
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(7*2.2) 


we  know  that 
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Y W1.,,  * dNN  h CR 

and  so  our  strategy  is  to  find  the  value  of  R at  which 
^dNN  YN  sufficien-tly  small. 

Of  course  €R  is  not  a known  quantity  since  it  in 
volves  the  exact  solution,  y(x).  Therefore  we  must  be  able 
to  approximate  €R  using  the  computed  solution,  Y. 

We  know  that 

€R  = a2Ny(xN+l)  * 02Ny  (xN+l* 

= ^a2N  4 ^2N^qN+l  " (7.2.4) 

The  value  of  yf*^^  must  be  approximated  in  (7.2.4).  We 
cannot  use  Y^+^  because  it  has  been  set  to  zero.  The  best 
value  available  to  us  is  YN.  We  approximate  ^ by  A.  Thus 
in  our  adaptive  algorithm  we  monitor  the  size  of 

6R  * *dNN  YN  ^a2N  4 02N(qN+l  “ A)}l. 

To  insure  that  6R  is  small  enough  so  that  no  sig- 
nificant error  is  introduced  we  require  that 

‘r  - ldjf  Yj  T(Xj) 1/100. 

This  factor  of  1/100  is  arbitrarily  chosen  and  can  be  con- 
sidered an  area  for  possible  tuning  of  the  algorithm. 

If  6r  is  not  small  enough,  a new  point,  xN+2#  is 
added  to  the  mesh  at 

xN+2  3 XN»1  4 h°o. 
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We  let  h,  start  with  the  value  of  h calculated  in 
the  first  iteration  of  the  algorithm.  Each  time  h^,  is 
used  after  that  its  size  is  doubled.  This  ensures  that 
R gets  large  quickly  enough  so  that  extra  points  are  not 
introduced  and  extra  iterations  not  performed.  This  factor 
of  2 is  again  an  area  for  possible  tuning  of  the  algorithm. 

The  algorithm  we  use  for  singular  problems  of  the 
first  type  is  as  follows « 

ALGORITHM  7.1 t This  algorithm  is  exactly  the  same 
as  Algorithm  5*6  but  we  include  the  following  statement 
in  Step  5i) . 

If  ,dNN  YN  ta2N  + 32N(Vu  ' A)  II  > EPS/100  then, 
a.  add  a point  xN+2  = xN+1  * h^. 

b*  h«  - 2h«» 

Results  for  problems  A - D of  Table  7.1  are  exam- 
ples of  Algorithm  7.1  in  operation. 

7.3  Singularity  in  q(x) 

When  the  function  q possesses  a singularity  at  a or 
b,  we  encounter  a similar  problem.  We  assume  for  most  of 
this  discussion  that  a is  the  singular  point  although  all 
arguments  are  valid  for  b as  well. 

It  is  obvious  that  the  point  a itself  cannot  be  an 
element  of  the  mesh  because  the  difference  equation  (4.3.1) 
would  require  q(a)  to  be  evaluated.  Therefore  we  approx- 
imate the  problem  by  letting  xQ  be  located  at  some  point 
close  to  a.  but  where  q(Xg)  is  defined. 
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Our  adaptive  algorithm  must  decide  if  xQ  is  close 
enough  to  a to  ensure  that  aX  is  small,  and  if  not  it  must 
introduce  a new  point  between  a and  xQ. 

The  error  caused  by  this  truncation  occurs  in  the 
first  equation  of  (7. 2.1).  We  again  have 

(A  - XB)y  = t(y)  - eQ 
where  now  eQ  • (€0,0, . . .,0)t  with 

€0  « any(x0)  ♦ B21yM(x0). 

As  in  Section  7.2,  we  approximate  by 
[a2i  + 021(qo  - and  we  monitor 

I d^^  Y1  £ ct 21  + ^21^0  ” A)}|. 

The  algorithm  we  use  to  accomplish  this  task  isi 

ALGORITHM  7.2»  This  algorithm  is  also  like  Algo- 
rithm 5*6  but  with  the  following  additions. 

1.  In  Step  1 the  basic  mesh  is  set  up  on  [x0,b]. 
The  user  can  supply  an  initial  value  for  xQ, 
otherwise  xQ  = a + (b-a)/8. 

2.  In  Step  5 include i 

If  ,dll  Ylta2l  + *21(q0  ’A)}!  > EPS/100 
then  split  xQ.  (add  point  at  (a+xQ)/2  ) 

Problems  C and  D of  Table  7.1  are  of  the  second 

type. 

The  inclusion  of  into  our  algorithm  does  not 


r 
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capture  all  of  the  error  introduced  by  the  truncation. 

Error  is  propagated  throughout  the  entire  eigenfunction  by 
this  truncation  process.  Techniques  for  estimating  this 
error  may  lead  to  even  better  error  estimates  than  those 
shown  in  Table  7.1. 

7.4  Numerical  Results 

In  this  section  we  show  results  of  running  Algo- 
rithms 7.1  and  7.2  on  several  singular  problems.  We  solve 
the  following  problems « 

Problem  Ai  (Morse  potential  problem) 
y + (X  - q(r))y  = 0 
y ( 0)  = 0,  y(r)  - 0 as  r * « 

where 

q(r)  = D[l  - exp(-a(r  - re))]2  - D 

with 

D = 188.4355.  a = 0.711248,  r = 1.9975. 

XQ  « -178.79850,  X4  « -110.80832. 

Problem  B> 

■ 2 

y (X  - x )y  » 0 

y(  - 30 ) ■ y(°°)  * 0 

XR  ■ 2n  ♦ 1,  n»0,l,2,... 

y (x)  « exp(-x2/2)H  (x) 
n n 

where  Hn(x)  is  the  Hermite  polynomial  defined  by«  HQ(x)  « 1, 
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H1(x)  • 2x,  3 2xHjt(x)  - 2kHjc_1(x),  k»0,l,... 

Problem  C« 

y ♦ U ♦ i)y  ■ o 

y(O)  = 0,  y(x)  - 0 as  x -•  *>. 

Xn  * -l/4(n+l) 2,  n-0,1, 2, . . . 

Problem  D» 

/ * U ♦ 1 - )y  - 0 

y(o)  = y(oo)  a 0 

* -l/4(L+n+l) 2,  n=0,l,2,... 

Table  7.1  shows,  for  each  of  the  above  problems, 
the  same  information  provided  by  Table  5*6.  In  addition 
we  present  the  calculated  values  of  a and  b. 

The  increased  accuracy  obtained  by  the  adaptive 
method  is  quite  evident  in  the  results  of  Table  7.1.  The 

values  of  a and  b,  used  to  compute  the  Numerov  results,  are 

those  calculated  by  the  adaptive  method.  Normally,  however, 
the  user  of  a uniform  mesh  method  must  select  his  own  a and 

b.  If  some  other  standard  values  are  chosen  for  the  prob- 

lems above  the  improved  accuracy  of  the  adaptive  method 
becomes  even  more  apparent. 


Table  7.1 
Singular  Problems 
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Mth  TOL  N 

a/b 

le| 

Problem  A,l0 

A 10“4  41 

0 

1.114D-4 

6.O83D-5 

5.0530-5 

N 41 

3075 

2.872D-4 

3.189D-4 

3.169D-5 

64 

4.058D-5 

5.436D-5 

1.378D-5 

Problem  A,X^ 

A 10‘4  270 

0 

1.918D-4 

3. 6950-5 

1 . 548D-4 

N 270 

14.625 

4.874D-3 

4.4950-3 

3.792D-4 

Problem  BtlQ 


A 10’6  136 

-8.75 

7.450D-7 

6.808D-7 

6.419D-8 

N 136 

8.75 

2.147D-6 

2.132D-6 

1.526D-8 

256 

1.707D-7 

1.704D-7 

3.4320-10 

Problem  B.Xg 

A 10"4  70 

-16.75 

1.1330-4 

8.6730-5 

2.658D-5 

N 70 

16.75 

1.084D-2 

9.222D-3 

1.618D-3 

128 

9.311D-4 

8.872D-4 

4.393D-5 

Problem  C#XQ 

A 10* 5 80 

1.097D-7 

2.702D-6 

2.530D-6 

1.7270-7 

N 80 

27.48125 

4.047D-3 

6.899D-4 

3.357D-3 

128 

1.706D-3 

3.094D-4 

1.396D-3 
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Table  7.1  (continued) 


Mth 

TOL 

N 

a/b 

let 

,ec> 

Problem  C, 

X2 

A 

10'5 

92 

1.097D-7 

177.10625 

5.513D-7 

5.327D-7 

1.849D-8 

N 

92 

2.252D-3 

2.032D-4 

2.049D-3 

128 

1.454D-3 

1.604D-4 

1.293D-3 

Problem  D, 

VL=3 

A 

i-* 

o 

i 

40 

.025 

336.7063 

4.114D-6 

3.326D-6 

7.880D-7 

N 

40 

9.935D-6 

1.184D-5 

1.910D-6 

64 

1.880D-6 

1.470D-5 

4.095D-7 

Problem  D, 

VL“3 

A 

10”5 

86 

.0125 

655.5906 

1.506D-6 

1.470D-6 

3.647D-8 

N 

86 

8.446D-6 

7.949D-6 

4.971D-7 

128 

1.803D-6 

1.492D-6 

3.110D-7 

The  diversity  of  the 

calculated 

a*s  and  b* 

s shows 

that  it  would  be  difficult  to  select  these  numbers  before- 


hand. A wide  range  of  values  for  h^  is  also  characteris- 
tic of  these  problems.  For  example,  in  Problem  C,  h^  takes 
on  values  from  1.907x10”?  to  4.9875  when  computing  *2. 
Similar  ranges  are  found  in  all  the  problems  of  this  chap- 
ter. No  user  can  be  expected  to  choose  such  a mesh  in 
advance . 


ft 


In  Figures  7.1  - 7.4  we  show  the  point  placement 
for  Problems  B and  C.  In  these  graphs  we  note  the  diver- 
sity of  the  h^'s  as  well  as  the  dense  mesh  in  steep  sec- 
tions of  the  functions  and  sparse  mesh  in  level  sections. 
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7.5  Order  of  Convergence 

When  the  problem  to  be  solved  has  a singularity  the 
order  of  convergence  results  of  Section  4.3  are  no  longer 
necessarily  valid. 

Let  us  consider,  for  example,  Problem  C.  The  adap- 
tive method  solves  the  problem  by  truncating  the  interval 
(0,«o)  to  (Xq.R).  We  now  solve  a perturbed  problem  which 
is  regular  on  (xQ,R).  We  would  obtain  O(h^)  convergence 
of  the  matrix  problem  to  the  perturbed  problem  as  h - 0, 
however,  convergence  of  the  perturbed  problem  to  the  orig- 
inal singular  S-L  problem,  as  h ■*  0,  xQ  -•  0,  and  R - «>, 
must  be  analyzed  separately  for  each  problem. 

The  convergence  rate  will  now  depend  on  hQ  = (xQ  - a) 
and  on  R.  In  this  example  we  find  that  while  solving  for 
we  have  |e|g  « 2.107X10"1  and  |eclg  * 2.095X10"1  with 
hQ  ■ .05.  After  convergence  we  have  |e|g0  = 2.702x10"^ 
and  |ecl80  = 1.727x10"?  with  hQ  = 1.907x10"?.  Computing 
a and  3 as  we  did  in  Section  5*3.4,  we  find  that  a » .903 
and  3 w 1.123.  These  figures  imply  that 

( A - A)  a 0(h) 
and  ( A - £)  » 0(h) • 

In  general,  these  values  of  a and  3 will  be  differ- 
ent for  each  singular  problem. 
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7.6  Calculation  of  a Transformation  Function 

Many  methods  for  solving  the  singular  S-L  problem 
involve  a mapping  of  the  problem  onto  some  other  domain. 

The  independent  variable,  x,  is  mapped  to  a new  variable 
r by  the  transformation 

r = f(x) . 

Examples  of  such  transformations  include  exponential 
• k x 

functions  (r  = e ),  log  functions  (r  = log  x),  and  others 
(e .g.  r • 

The  purpose  of  such  a transformation  is  to  alter  the 
differential  equation  so  that  a uniform  mesh  method  can  be 
used  effectively.  As  an  example  we  look  at  the  hydrogen 
atom  equation,  Problem  C of  Chapter  ?.  We  know  that  a so- 
lution can  only  be  found  efficiently  if  the  mesh  is  fine 
near  zero  and  very  coarse  for  larger  x.  The  transforma- 
tion f(x)  should  thus  spread  out  the  smaller  values  while 
compressing  larger  ones.  The  adaptive  method  of  this 
thesis  is  useful  in  selecting  such  a transformation. 

We  want  a function  f(x)  with  the  property  that 

t 

Ar  = f (x)ax  » constant. 

We  solve  the  example  problem  using  the  adaptive  algorithm. 
This  gives  us  an  optimal  mesh  of  ax's.  We  now  look  for  a 
function  which  satisfies 

• 

f (x)  Ax  as  constant. 
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For  this  problem  we  found 
f(x)  = jij 

to  be  an  appropriate  function.  Table  7.7  shows  values  of 
f (x)Ax  for  every  x at  which  a step  size  change  occurs. 

N was  computed  to  be  67  for  this  problem.  Thus  the 
transformed  problem  should  use  Ar  = ^ a 1.4925D-2.  We 
see  that  many  of  the  values  in  Table  7.2  are  close  to  this 
number.  Discrepancies  are  due  to  the  discontinuity  of  our 
computed  mesh. 

Such  analysis  would  need  to  be  done  for  each  prob- 
lem and  selection  of  a proper  f(x)  could  be  difficult.  The 
results  of  this  section  show,  however,  that  such  a use  of 
the  adaptive  method  can  be  worthwhile  for  selecting  a 
transformation  if  a uniform  mesh  method  must  be  used. 


Table  7.2 

f (x)ax  for  Problem  C 


X 

f*  (x)ax 

1.250D-2 

1.219D-2 

5.000D-2 

3. 534D-2 

1.279D-1 

1.531D-2 

2.838D-1 

1.494D-2 

6.734D-1 

5.566D-2 

1.609D0 

4.581D-2 

5. 66 IDO 

1.405D-2 

1.626D1 

4.186D-3 

2.249D1 

4.518D-3 

1 


Chapter  8 


PERFORMANCE  EVALUATION 

In  this  chapter  we  evaluate  the  performance  of  the 
adaptive  algorithm  in  terms  of  space  and  time  utilization. 
We  show  relationships  between  the  index  of  the  desired  eig- 
envalue and  the  space  and  time  requirements. 

8.1  Space  Requirements 

The  measure  of  space  requirements  of  our  adaptive 
algorithm  is  N,  the  number  of  points  in  the  final  mesh. 

We  know  from  the  many  problems  of  Chapters  5 and  7 
that  for  a given  problem  N is  a function  of  both  i,  the 
index  of  the  desired  eigenvalue,  and  TOL,  the  allowable 
error  in  the  numerical  solution.  In  order  to  simplify  the 
analysis  we  assume  that  TOL  is  constant  and  we  exhibit 
the  behavior  of  N as  i increases.  Thus  we  define  the 
function  S(i)  = N. 

Figures  8.1  - 8.4  show  the  computed  value  of  N plot- 
ted against  i for  several  different  problems.  TOL  is  held 
constant  at  10  . 

8.1.1  Regular  Problems 

The  function  S is  different  for  each  problem  and  so 
a separate  analysis  must  be  performed  on  each  to  obtain  a 
closed  expression  for  S(i).  Figures  8.1  - 8.3  suggest, 
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however,  that  there  is  a pattern  to  the  behavior  of  reg- 
ular problems. 

We  look  at  the  simplest  of  these  problems  to  analyze, 
Problem  I.  We  know  from  Chapter  4,  (4.3.22),  that  when  our 
adaptive  method  is  used  to  solve  this  problem  for  we  have 

» - Y^D_1D_1  t (y)  / YtD"1D“1BY. 

This  is  the  quantity  we  want  to  be  constant  for  all  i. 

We  let  G = 1 / YtD’1D'1BY  and  recall  from  Chapter 
4 that  G = 0(h”^).  So  we  have 


» - GYtD"1D'1  C(y) 

N+l  _ 

= - G I Y.d:?  t(x,). 


0=0 


■ns 


Now,  because  we  have  assumed  that  the  mesh  is  uni- 
form almost  everywhere,  we  know  that 


t(x.)  = 0(h6)y(vi) (x.)  + 0(h8) . 

J J 

We  recall  Problem  I to  be 

y + X y = 0 (8.1.1) 

y(0)  = y ( 1 ) = 0 

k^  = (i+l)2tr2  , i=0,l,2,... 


Differentiating  (8.1.1)  4 times  we  get 
y*vi^(x)  = - k ^y(x) 


and  so  we  have 
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a Q 


N+l  -2 

£0  Vjj  0(h°)X^y(Xj) 


,<u3. 


Since  we  know  the  eigenvalues  of  the  problem,  we 

3 6 6 6 

can  replace  X£  by  (i+1)  it  or  simply  0(i  ).  Thus 

N+l  ? a * 

AXi  s G L Yjd“J  0(h°)  0(i°)  y(Xj).  (8.1.2) 


But  (8.1.2)  is  a sum  of  N ♦ 2 or  0(N)  terms  so  we  have 


AX i « G 0(N)  0(h6)  0(i6). 

Now  replacing  G by  0(h_1)  and  O(h^)  by  0(N"^)  (because 
h ■ 0(l/N)}  we  have 


AXi  a 0(i6)  0(N‘4). 


This  says  that  for  AX^  to  be  constant  for  all  i,  N 
">/2 

must  increase  like  LJ/  . The  solid  lines  on  Figures  8.1 
8.3  show  functions  of  the  form 


N = S ( i ) = c0  + Cl  i3/2  (8.1.3) 

fitted  to  the  data  by  least  squares. 

Actually  each  regular  problem  must  be  analyzed 
individually  to  get  such  a relationship*  However,  we  note 
that  since  eigenfunctions  of  all  regular  problems  are  sim- 
ilar in  shape  we  might  expect  the  S(i)  to  be  similar  also. 
This  hypothesis  is  given  even  more  credence  by  observing 
the  excellent  fits  of  the  least  squares  functions  in  Fig- 
ures 8.2  and  8.3. 


The  exponent,  3/2,  in  equation  (8.1. 3)  is  a result 
of  several  approximations  or  assumptions.  Because  of  this, 
we  cannot  place  too  much  faith  in  the  exact  value  of  3/2. 

We  have,  however,  fitted  linear  and  quadratic  least  squares 
functions  to  the  data  of  Figure  8.1.  When  we  compare  the 


values 


L (N  - S(i)  )' 
i -0 


for  the  3 different  functions,  S(i),  we  find  that  the 
value  for  (8.1.3)  is  smallest.  ( 3709.0  compared  to  3710.9 
and  5087. 9).  Thus,  while  3/2  may  not  be  the  exact  exponent, 
it  does  appear  that  the  value  lies  between  1 and  2. 

8.1.2  Singular  Problems 

The  process  used  above  was  also  performed  on  a sing- 
ular problem.  Problem  B. 


y ♦ U - x2)y  = 0 

y(-  30 ) - y(«>)  = 0 

Xn  = 2n  + ! , n-0 ,1,2, ... 


The  problem  was  solved  for  i=0,l,...,6  using  Algorithm 
7.1,  holding  TOL  constant  at  10-2.  Figure  8.4  shows  the 
plotted  results.  The  data  appears  to  be  linear  with  re- 
spect to  i.  The  solid  line  is  a least  squares  straight 
line.  The  computed  values  of  a and  b range  from  (-3.5,3. 5) 
to  (-4.75,4.75).  Each  was  started  at  (-1.0, 1.0). 


158 


8.1.3  Conclusions  of  Space  Requirement  Findings 

We  found,  for  Problem  I and  possibly  many  other  reg- 
ular problems,  that  while  space  does  increase  exponentially 
with  respect  to  i,  it  is  not  at  a prohibitive  rate. 

Using  the  least  squares  function  N = S(i)  = 18.136  + 
7.921  i^2  we  find,  for  example,  that  A.?0  of  Problem  I 
could  be  computed  with  approximately  727  mesh  points  to  an 
absolute  accuracy  of  10  . Since  ~ **352.5  this  repre- 

sents a relative  accuracy  of  about  2.3*10"^.  In  the  sing- 
ular problem  we  observed  a linear  relationship  between  N 

and  i.  In  that  problem  could  be  computed  with  about 

327  points  to  a relative  accuracy  of  2.4xl0-i|. 

We  note,  however,  that  the  number  of  points  required 

for  a singular  problem  changes  when  the  initial  a and  b 
are  changed.  If  the  final  a and  b were  known  before  run- 
ning the  algorithm  they  could  be  input  thus  reducing  the 
size  of  the  mesh. 

8 . 2 Time  Requirements 

The  second  measure  of  the  performance  of  our  method 
is  time.  Figures  8.5  - 8.8  show  the  actual  CPU  time  re- 
quired to  compute  each  eigenvalue.  Our  test  program  is 
written  in  FORTRAN  and  compiled  with  the  FORTRAN  G compil- 
er on  an  IBM  370/168.  The  time  was  measured  by  calling 
system  clock  routines  at  the  start  and  at  the  end  of  the 
computation.  The  times  shown  are  only  for  the  portion  of 
the  program  corresponding  to  the  adaptive  algorithm.  Ini- 
tialization and  input/output  are  excluded.  All  times  are 
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in  milliseconds.  Exact  values  of  the  times  are  not  what  we 
wish  to  show  but  rather  the  relative  times  as  i increases. 

All  the  processes  done  in  the  adaptive  algorithm 
are  known  to  be  linear  in  time  with  respect  to  N.  This 
includes  matrix  reduction,  mesh  refinement,  and  coeffi- 
cient and  eigenvalue  computation.  Thus  if  times  were 
measured  for  only  the  last  iteration  of  each  data  set  we 
would  see  N and  T (time)  related  linearly  while  i and  T 

*1 /p 

would  again  be  related  like  T = ).  However,  for  the 

algorithm  to  arrive  at  the  higher  values  of  N it  must  com- 
pute eigenvalues  for  several  lower  values  of  N first.  Each 
of  these  computations  in  turn  takes  several  Rayleigh  quo- 
tients iterations.  Thus  as  i increases,  thereby  increasing 
N,  the  time  shown  is  a sum  of  several  iterations. 

To  analyze  this  relationship  we  define  a new  var- 
iable M to  be  the  sum  of  N x (number  of  Rayleigh  quotient 
iterations)  summed  over  each  N used  to  compute  X.  For 
example,  in  computing  X^  of  Problem  I we  find  that  matrix 
problems  of  sizes  N = 8,  16,  32,  38,  50,  52,  66,  78,  84, 
and  139  must  be  solved  and  that  respectively  4,  4,  3,  4,  4, 
4,  3«  3#  and  3 Rayleigh  quotient  iterations  are  re- 

quired. Thus  M = 1887. 

If  the  process  is  indeed  linear  with  respect  to 
time  within  each  iteration  then  we  must  expect  that  M and 

v 

T are  linearly  related.  Figure  8.9  shows  M and  T plotted 
for  Problem  I.  The  solid  line  is  the  least  squares  line 
T = 1.580  M ♦ 73*9^0.  Similar  linear  relationships  exist 
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for  the  other  example  problems  of  this  chapter. 

All  the  data  in  Figures  8.5  - 8.9  were  obtained 
by  initially  letting  N = 8.  This  value  is,  however,  an 
input  to  the  algorithm.  If  the  user  wants  to  decrease  the 
time  for  higher  eigenvalues  he  can  set  N = 16,  32,... 
or  some  higher  value  based  on  his  knowledge  of  the  prob- 
lem, This  would  eliminate  some  of  the  earlier  iterations 
thus  reducing  the  time  required  for  convergence.  We  must 
note,  however,  that  these  early  iterations  use  much  less 
time  than  the  later  ones  because  N is  still  small.  Thus 
the  time  savings  would  possibly  not  be  a considerable  one. 
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N vs  i - Problem  B 
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Figure  8.5 
T vs  i - Problem  I 
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Figure  8.8 
T vs  i - Problem  B 
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T vs  M - Problem  I 
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8.3  Comparison  With  an  Initial  Value  Method 

In  this  section  we  present  some  numerical  results 
comparing  our  algorithm  with  SLEIGN,  the  software  pack- 
age described  in  [l]  and  briefly  discussed  in  Section  3*5 
of  this  thesis. 

We  compared  performance  on  regular  problems  by  run- 
ning both  algorithms  on  several  sample  problems.  Compar- 
ative times  (T),  measured  in  milliseconds,  are  shown  in 
Table  8.1  along  with  actual  errors  (|e|)  and  error  esti- 


mates (|eegtl).  Another  important  measure  of  performance 
is  the  value  1®est^e^*  ra-ti°  of  error  estimate  to 
actual  error.  Ideally  this  ratio  should  be  equal  to  1. 

Any  value  significantly  less  than  1 denotes  a small  esti- 
mate which  would  have  caused  the  algorithms  to  stop  pre- 
maturely. A value  greater  than  1 indicates  that  the  algo- 
rithm did  extra  work  by  making  the  actual  error  too  small. 

Data  from  SLEIGN  for  comparing  performance  on  sing- 
ular problems  was  taken  directly  from  [l].  Times  are  not 
available  so  we  look  only  at  Problems  B,  C, 

and  D were  run  by  each  algorithm  and  Table  8.2  shows  the 
comparison.  In  addition,  Table  8.3  shows  values  for  all 
other  eigenvalues  computed. 

TOL  was  set  at  10*^  for  problem  I.  For  all  other 

-4 

problems  TOL  was  set  at  10 


Table  8.1 
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Adaptive  vs  SLEIGN  - Regular  Problems 


Prob  Meth 

N 

T 

let 

*eest* 

Ratio 

'•o' 

i.*.0 

A 

44 

173^ 

2.395D-6 

2.393D-6 

.9992 

1.732D-9 

S 

- 

1767 

3.052D-5 

9.998D-6 

.3276 

- 

II. *0 

A 

16 

150 

1.431D-5 

1.350D-5 

.9434 

8.127D-7 

S 

- 

844 

2.861D-6 

6.337D-6 

2.215 

- 

n,x2 

A 

68 

894 

3.129D-5 

3.114D-5 

.9952 

1.533D-7 

S 

- 

5486 

1. 526D-5 

2.827D-5 

1.853 

- 

III.‘o 

A 

20 

306 

4.06 2D- 5 

4.051D-5 

.9973 

1.179D-7 

S 

• 

1274 

7 .40 3D- 5 

2.520D-5 

.3404 

• 

Table 

8.2 

Adaptive  vs  SLEIGN  - Singular  Problems 
Computed  by  Both  Algorithms 

Problem 

|eestl/lel  “ Adaptive 

'Set'/'o 

| - SLEIGN 

.9138 

.3100 

C.A0 

.9396 

10 

.00 

d.a0 

.8085 

.2375 

I 


Table  8.3 

Adaptive  vs  SLEIGN  - All  Singular  Problems 
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1. 


1 4 

Prob 

Adaptive 

Ratio(XQ)  Ratio(A2) 

Prob 

SIEIGN 

Ratio(XQ) 

Ratio(X^Q) 

u 

A 

.5^61 

.1926 

1 

.3100 

.7988 

r-. 

B 

.9138 

.7655 

2 

.2375 

5.02 

C 

.9363 

.9663 

3 

.1850 

mm 

D 

.8085 

.9761 

4 

162.79 

69.09 

u 

5 

6.895D-3 

1.063D-2 

6 

1.200 

186.49 

13 

10.00 

1.135 

(Note  1 Values  are  shown  in  Table  8.3  only  for  those  prob- 
lems for  which  the  actual  error  is  known  to  be  reliable.) 

This  limited  comparison  appears  to  show  that  the 
adaptive  method  is  preferable  to  this  initial  value  method 
for  two  reasons.  First,  the  error  estimate  obtained  by  our 
method  is  much  more  reliable  than  the  one  of  the  initial 
value  method.  Values  of  !eest1/lel  are  consistantly  around 
1.0  for  our  method,  ranging  from  .1926  to  .9992,  while  the 
ratios  for  SLEIGN  vary  from  6.895*10”^  to  162,791,  Second, 
the  method  appears  to  be  significantly  faster  for  similar 
accuracies.  It  must  be  noted,  however,  that  SLEIGN  solves 
the  problem  without  requiring  it  to  be  put  in  the  normal 
form. 
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Chapter  9 


SUMMARY  AND  CONCLUSIONS 

In  this  thesis  we  developed  an  adaptive  finite- 
difference  method  for  solving  numerically  the  Sturm- 
Liouville  problem.  The  Sturm-Liouville  theory  is  pre- 
sented in  Chapter  2.  Chapter  3 serves  as  background  on 
the  many  other  techniques  available  for  solving  this  prob- 
lem. We  found  that  most  current  work  is  done  with  some 
variation  of  the  initial  value  or  shooting  method.  These 
methods,  however,  traditionally  use  Newton's  method  to 
find  a root  or  eigenvalue.  This  leads  to  quadratic  con- 
vergence to  the  desired  eigenvalue. 

In  Chapter  4 we  presented  a non-uniform  finite- 
difference  method.  We  showed  that  when  h,  the  maximum 
step  size,  is  sufficiently  small  this  method  is  capable 
of  computing  all  eigenvalues  and  that  cubic  convergence 
is  experienced  via  Rayleigh  quotient  iteration.  We  showed 
also  that  with  one  application  of  the  deferred  correction 
algorithm,  the  method  is  O(h^)  in  the  eigenvalue  and  eig- 
envector. 

Two  algorithms  were  presented  in  Chapter  5.  The 
first  is  a non-adaptive , iterative  algorithm.  It  shows 
the  capability  of  a non-adaptive  method  which  uses  a uni- 
form mesh. 
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The  second  algorithm  is  our  adaptive,  iterative 
method.  It  uses  the  non-uniform  finite-difference  method 
of  Chapter  4.  The  D“^  Solution-Weighted  strategy  is  found 
to  be  the  best  for  our  problem  and  it  is  used  in  the  algo- 
rithm. We  found  that  the  adaptive,  iterative  algorithm 
does  indeed  obtain  accuracy  equivalent  to  the  uniform  mesh 
method  while  using  fewer  grid  points. 

In  Chapter  6 we  solve  the  problem  of  finding  a 
starting  value  for  the  method  and  for  monitoring  the  re- 
cycled eigenvalue  from  iteration  to  iteration.  We  present- 
ed an  efficient  method  for  estimating  the  desired  eigen- 
value close  enough  to  ensure  proper  convergence.  Sturm 
sequences  are  used  to  accomplish  these  tasks.  Thus,  our 
method  is  self-starting,  making  it  more  practical  to  use 
than  most  other  current  methods. 

Singular  problems  were  addressed  in  Chapter  ?•  We 
showed  that  with  a few  additions  to  the  adaptive  algo- 
rithm, singular  problems  are  solved  quite  well.  We  found 
that  an  adaptive  method  is  essential  for  the  efficient  sol- 
ution of  these  problems. 

Finally,  in  Chapter  8,  we  showed  the  relationship 
between  the  index  of  the  eigenvalue  and  the  space  and  time 
requirements  of  our  algorithm.  We  found  that  our  method  is 
practical  for  use  in  computing  more  than  just  the  lowest 
eigenvalues.  We  noted,  however,  that  these  relationships 
are  highly  problem  dependent.  We  also  showed  that  our  meth- 
od can  outperform  an  existing  initial  method  in  terms  of 


accuracy  of  error  estimates  and  speed  of  computation. 

We  conclude  that  the  automatic,  adaptive  Sturm- 
Liouville  solution  method  developed  in  this  thesis  is  an 
essential  tool  for  solving  the  Sturm-Liouville  problem. 

It  makes  the  finite-difference  technique  competitive  with, 
or  preferable  to,  the  other  methods  being  used. 
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